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Observational advances over the last decade have enabled high-resolution, interferometric 
studies of forming multiple systems, statistical surveys of multiplicity in star-forming regions, and 
new insights into disk evolution and planetary architectures in these systems. In this review, we 
compile the results of observational and theoretical studies of stellar multiplicity. We summarize 
the population statistics spanning system evolution from the protostellar phase through the 
main-sequence phase and evaluate the influence of the local environment. In short, most stars are 
born in multiple stellar systems, most main sequence stars are members of multiple systems, but 
most star systems are single. We describe current models for the origin of stellar multiplicity and 
review the landscape of numerical simulations and evaluate their consistency with observations. 
We review the properties of disks and discuss the impact of multiplicity on planet formation 


and system architectures. Finally, we summarize open questions and discuss the technical 
requirements for future observational and theoretical progress. 


1. INTRODUCTION 


The formation of multiple star systems — systems of 
two or more gravitationally bound stars with separations 
S 0.1 pc - takes place during the earliest phases of star for- 
mation. The majority of such systems form and evolve to 
their final configuration during the time period spanned by 
the collapse of dense cores through the end of mass accre- 
tion. This review summarizes the current statistics of mul- 
tiple stars, describes theoretical models and observational 
evidence for different modes of multiple star formation, and 
discusses the implications for planet formation and system 
evolution. 

Studying multiple formation has historically been chal- 
lenging from several perspectives. Dust in star-forming 
regions absorbs short-wavelengths of light and re-emits at 
longer wavelengths, requiring high-resolution, radio obser- 
vations to detect forming stars. Theoretical models, mean- 
while, require both large dynamic range (^1 au -10 pc) 
and treatment of non-linear, complex physical processes, 
including turbulence, radiative feedback, and magnetic ef- 
fects. Consequently, when the last Protostars and Planets 
conference occurred in 2013, the mechanisms for multi- 
ple formation were largely theoretical and the properties of 
young multiple systems were poorly constrained. 

At that time a major expansion of the Very Large Array 
(VLA) had just concluded, and only the first cycle of ob- 
servations by the Atacama Large Millimeter/submillimeter 
array (ALMA) had been taken. Since then, the power of 
these telescopes has revolutionized our understanding of 
young forming multiples, and stunningly detailed images of 
disks have provided a comprehensive understanding of the 
impact of binaries on disks, both in a statistical sense and 
within individual multiple systems. Meanwhile large-scale, 
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uniform surveys of star-forming regions, clusters, and field 
stars have generated new insights into the statistics of stellar 
multiplicity and its dependence on mass, age, metallicity, 
and cluster density. For already-formed planetary systems, 
radial velocity (RV) and high-resolution imaging surveys of 
Kepler and TESS planet hosts have recently revealed that bi- 
naries substantially sculpt planet formation, shedding new 
light on planet formation and migration models. Compu- 
tational advances in speed and methodology have enabled 
studies of multiple star formation that span an unprece- 
dented dynamic range, model different star cluster environ- 
ments, and explore the role of multiple physical effects act- 
ing in concert. There is now consensus that most stars both 
form and reside on the main sequence with at least one stel- 
lar companion: single stars are the exception not the rule. 

This review begins by introducing key definitions and 
summarizing the statistics of observed multiple star sys- 
tems ($2), with a special focus on their configuration and 
properties at the end of the star formation process. We then 
proceed with an overview of theoretical models describing 
the origin and early evolution of multiple star systems in 
$3. Using the framework of these models, in $4 we exam- 
ine protostellar multiple systems, assessing how well they 
confirm or refute theoretical expectations. In $5 we discuss 
the impact of stellar multiplicity on planet formation and 
system architectures. $6 summarizes remaining open ques- 
tions and discusses observational and theoretical prospects 
for enhancing the understanding of multiple systems over 
the coming decade. We conclude in $7. 


2. OBSERVED STELLAR MULTIPLICITY 


Observations show that stellar multiplicity evolves sub- 
stantially during formation and changes more slowly there- 
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after. In order to distinguish between primordial multiplic- 
ity, the results of dynamical evolution, and the final system 
configurations for main-sequence (MS) field stars, we sep- 
arately consider the multiplicity of very young stars, i.e., 
protostars, and older sources. We adopt the standard con- 
vention of grouping young stellar objects (YSOs) according 
to their spectral energy distribution, i.e., by “Class” as de- 
fined by Lada and Wilking (1984) and Adams et al. (1987), 
whereby younger, more embedded sources have more red- 
dening from an envelope and disk. We caution, however, 
that the inclination of the source outflow and disk with re- 
spect to the line of sight affect the degree of extinction and 
may shift an older source from a later to an earlier Class 
or vice versa (Robitaille et al. 2007; Offner et al. 2012a). 
Thus, while Class provides a rough, if imperfect, measure 
of sample age it is an unreliable age indicator for individual 
sources. 

In this section, we compile and compare the multiplicity 
statistics of MS and non-embedded pre-MS stars (mostly 
Class II/III T Tauri stars) from several surveys (mostly op- 
tical and near-IR). These statistics represent the final out- 
come of the star formation process, providing insight into 
how stars and planets form and must subsequently evolve, 
and serve as a critical benchmark for simulations and mod- 
els. A detailed analysis of embedded Class 0/I binaries ob- 
served by the VLA and ALMA is discussed in $4.2. 

Many multiplicity properties also vary as a function of 
the physical distance between sources. In the literature, pair 
separations are often characterized as "close" or “wide,” 
however, these terms have no fixed definition. Throughout 
this review, we define close, intermediate, and wide sepa- 
rations as < 10 au, ~ 10 - 300 au and = 300 au, respectively, 
unless otherwise specified. 


2.1. 


No single observational technique provides a complete 
picture of multiplicity. For MS and T Tauri stars, spectro- 
scopic and eclipsing binaries probe close separations within 
a < 10 au (Abt et al. 1990; Latham et al. 2002; Melo 
2003; Sana et al. 2012; Moe et al. 2019; Kounkel et al. 
2019). Interferometry, sparse aperture masking, Hubble 
Space Telescope (HST) imaging, and astrometry reveal bi- 
naries across close to intermediate separations of a = 1 - 
300 au (Reid et al. 2001; Gizis 2002; Kraus et al. 2008; 
Raghavan et al. 2010; Dieterich et al. 2012; Rizzuto et al. 
2013; Sana et al. 2014; De Furio et al. 2019). Speckle 
imaging, lucky imaging, and adaptive optics (AO) resolve 
binaries across intermediate to wide separations of a = 10- 
1,000 au (Shatsky and Tokovinin 2002; Close et al. 2003; 
Law et al. 2008; Bergfors et al. 2010; Janson et al. 2012; 
De Rosa et al. 2014; Ward-Duong et al. 2015; Tokovinin 
and Bricefio 2020). Common proper motion confirms that 
wide binaries with a — 300 - 20,000 au are in fact gravita- 
tionally bound (Lépine and Bongiorno 2007; El-Badry and 
Rix 2018; Winters et al. 2019; Hartman and Lépine 2020). 
Other techniques can identify unresolved binaries, e.g., ex- 
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cess luminosity method or deblending of spectroscopic bi- 
naries, but they cannot determine their orbital separations 
(Hubrig et al. 2001; Bardalez Gagliuffi et al. 2014; Gullik- 
son et al. 2016; Widmark et al. 2018). For volume-limited 
populations of M-dwarfs (Winters et al. 2019) and solar- 
type MS stars (Raghavan et al. 2010; Tokovinin 2014b), 
the combination of these overlapping techniques provides 
a >90% complete census of MS companions orbiting the 
primary. 

We first define the multiplicity fraction (sometimes 
called the binary fraction) as the fraction of primaries with 
at least one companion: 


Mi Bep Tere. 
S+B+T+Q+... 


where S, B, T, and Q are the number of single stars, bina- 
ries, triples, and quadruples, respectively. Higher ordered 
multiples are rare but exist, including the two known sep- 
tuples, AR Cas and v Sco, both of which have early-B 
primaries (Eggleton and Tokovinin 2008). Similarly, the 
triple/high-order fraction is: 


(1) 


T-Q-.. 


THF = 
S+B+T+Q+.. 


(2) 


Meanwhile, the companion frequency is the average fre- 
quency of companions per primary: 


B+2T +30 +... 


CF = 
S+B+T+Q+... 


= | donedissse 43 


where fiog is the frequency of companions per decade of 
orbital separation. 

Measurement uncertainties in the multiplicity fraction 
and triple/high-order fraction follow binomial statistics, 
while Poisson statistics govern the companion frequency, 
which can exceed unity (Burgasser et al. 2003; Sana et al. 
2014; Moe and Di Stefano 2017). Here we report multiplic- 
ity statistics with respect to volume-limited samples. For 
magnitude-limited samples, corrections for Malmquist bias, 
which is sometimes called the Branch bias in the context of 
binary stars (Branch 1976), must be applied since unre- 
solved twin binaries are Am = 0.75 mag brighter than their 
single-star counterparts. Systematic uncertainties, such as 
those associated with correcting for incompleteness, gener- 
ally dominate the overall error for samples exceeding a few 
hundred systems. In the following subsections, we compute 
bias-corrected ME, THF, and CF from various surveys and 
present the results in Table 1 and Fig. 1. 

We display the separation distribution fioga of all com- 
panions for different stellar populations in Fig. 2. To dis- 
tinguish the effect of hierarchical triples on the overall sep- 
aration distribution, we define à4j and a;,, as the median 
separations of all companions and of inner binaries only 
(excluding outer tertiaries), respectively. Close binaries ex- 
hibit different trends with respect to primary mass, mass ra- 
tio, metallicity, and environment compared to intermediate 
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Fig. 1.—: Bias-corrected multiplicity fraction (left; thick), triple/high-order fraction (left; thin), and companion frequency (right) of 


BDs and MS stars. All three quantities increase monotonically with primary mass. The indicated spectral types at the top roughly 
correspond to the mean primary masses of field dwarfs. See Table 1 for detailed references. 


and wide separation binaries (see $2.2.7). We therefore de- 
fine the close binary fraction CBF and wide binary fraction 
WBF as the fraction of primaries with at least one com- 
panion below a « 10 au and across a = 100- 10,000 au, 
respectively. We display both Gay and à; in Fig. 3 and both 
CBF and WBF in Fig. 4 as a function of primary mass Mı 
and other stellar parameters. 

The minimum companion mass Mcomp or mass ratio 
q = Meomp/My included in reported multiplicity statistics 
varies from survey to survey, but typically brown dwarf 
(BD) primaries include all BD companions, FGKM-dwarfs 
include only MS companions above Mcomp > 0.075 Mo, 
and OBA stars include MS companions above q > 0.1. 
About 2096 of companions within a « 10 au of AFG pri- 
maries are white dwarfs (WDs; Moe and Di Stefano 2017; 
Murphy et al. 2018), and nearly 30% of OB stars are the 
products of binary evolution (de Mink et al. 2013). In order 
to compare field MS multiplicity statistics to their pre-MS 
counterparts, we exclude evolved stars or systems with WD 
companions when possible. 

Historically, the mass-ratio distribution has been approx- 
imated by a single power-law f, c q7 across 0 < q < 1, 
but larger samples indicate that up to three parameters may 
be required to accurately fit the distributions (Duchéne and 
Kraus 2013; Moe and Di Stefano 2017; El-Badry et al. 
2019). For this review, we wish to easily compare trends, 
even if the samples are small or incomplete toward small 
mass ratios. We therefore fit a simplified power-law slope 
Ytrunc across the truncated interval 0.4 < q < 1 where all 
the compiled surveys are complete (see Table 1). 


2.2. Multiplicity Statistics versus Primary Mass 
2.2.1. FGK dwarfs 


Frequencies: About half of field FGK primaries have 


MS companions, which follow a broad lognormal separa- 
tion distribution peaking near à4j = 40 au (Duquennoy and 
Mayor 1991; Raghavan et al. 2010; Tokovinin 2014a, see 
Figs. 1-4). In addition to plotting the separation distribu- 
tion directly from the Raghavan et al. (2010) data in Fig. 2, 
we also display a lognormal fit to fioga, where the fit pa- 
rameters are listed in Table 2. 

Across close and intermediate separations, the binary 
fraction of FGK MS stars in young open clusters are consis- 
tent with their field MS counterparts (Bouvier et al. 1997; 
Patience et al. 2002; Elliott et al. 2014; Deacon and Kraus 
2020; Torres et al. 2021). The CBFs measured in old open 
clusters also match the field values, albeit the cluster cores 
exhibit a substantial excess of binaries due to mass segrega- 
tion (Geller and Mathieu 2012; Geller et al. 2021b). Mean- 
while, Gaia observations reveal a deficit of wide binaries in 
open clusters, especially those with higher stellar densities, 
which is likely due to dynamical disruptions (Deacon and 
Kraus 2020; Niu et al. 2020). 

Mass Ratios: The overall mass-ratio distribution of 
solar-type binaries is roughly uniform, but closer solar-type 
binaries systematically favor larger mass ratios (Duquen- 
noy and Mayor 1991; Raghavan et al. 2010; Tokovinin 
2014b; Moe and Di Stefano 2017). Specifically, FGK bi- 
naries across close and intermediate separations exhibit a 
30% and 10% excess fraction of twins with q > 0.95, re- 
spectively, relative to the underlying uniform distribution 
(Lucy and Ricco 1979; Tokovinin 2000; Moe and Di Stefano 
2017). These observations suggest that solar-type binaries 
within a « 200 au coevolved in a shared mass reservoir, 
perhaps a circumbinary disk (see $3.2). El-Badry et al. 
(2019) also found that wide binaries exhibit a very small 
but statistically significant excess twin fraction, suggesting 
some wide companions originate from intermediate separa- 
tions but then subsequently widen dynamically. Aside from 
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Fig. 2.—: Frequency fioga of companions per decade of orbital separation. Left: the separation distribution of MS binaries as a 
function of spectral type and metallicity. Late-M (red; Winters et al. 2019), early-M (magenta; Winters et al. 2019), and FGK (black; 
Raghavan et al. 2010) binaries follow lognormal distributions (dashed fits - see parameters in Table 2). The distribution of all companions 
to mid-B (green; Moe and Di Stefano 2017) and O (blue; Sana et al. 2012, 2014) stars roughly follow Opik’s law. The dotted blue line 
shows the fit to inner binary companions to O stars, which skew significantly toward shorter separations. The close binary fraction of 
FGK dwarfs increases significantly at low metallicity (orange; Moe et al. 2019). Right: separation distribution of young solar-type 
binaries across different star-forming environments and classes, where we display the FGK MS distribution (black) for reference. The 
MS distribution is well matched by the bias-corrected close binary fraction within a < 10au of Class II/III T Tauri stars inferred from 
spectroscopy (green; Kounkel et al. 2019). However, imaging surveys reveal an excess of companions to pre-MS stars across a = 10- 
50 au in most environments, in conflict with the MS and Kounkel et al. (2019) results. Across a = 100 - 5,000 au, the dense Orion Nebula 
Cluster exhibits a deficit of wide companions (magenta; Reipurth et al. 2007; Duchéne et al. 2018), intermediate density regions like 
Upper Scorpius are consistent with the field population (blue; Tokovinin and Briceño 2020), and sparse, low-density Taurus shows a 
slight excess (orange; Kraus et al. 2011). Class I and flat-spectrum protostars (cyan) and especially Class 0 protostars (red) in Perseus 
and Orion exhibit a substantial excess of companions across a = 1,000 - 10,000 au (Tobin et al. 2022). Note that the two panels depict 


populations with very different ages, underscoring the evolution in multiplicity that occurs during formation. 


the excess twins, solar-type binaries beyond a > 200 au are 
skewed toward smaller mass ratios, but their distribution is 
still top-heavy compared to random pairings drawn from 
the IMF (Lépine and Bongiorno 2007; Moe and Di Stefano 
2017; El-Badry et al. 2019). 

There is a dearth of BD companions within a « 0.5 au 
to solar-type primaries, which is commonly called the BD 
desert (Grether and Lineweaver 2006; Csizmadia et al. 
2015; Shahaf and Mazeh 2019). Meanwhile, across inter- 
mediate separations, the roughly uniform mass-ratio distri- 
bution of solar-type binaries extends down to the BD-planet 
boundary near Mcomp = 13 Mj (Wagner et al. 2019; Nielsen 
et al. 2019). While the multiplicity fractions of MS stars 
reported in Table 1 exclude BD companions, only ~4% 
of solar-type stars have BD companions, affecting the in- 
tegrated multiplicity statistics very little. 

Eccentricities: Eccentricity distributions are typically 
referenced to a thermal distribution defined as p, = 2e, 
where pe is the probability of given eccentricity e (Am- 
bartsumian 1937; Heggie 1975). The average e of solar- 
type binaries increases with orbital separation, transitioning 
from a sub-thermal (smaller eccentricities) to thermal dis- 
tribution distribution near 200 au and to super-thermal be- 


yond » 1,000au (Tokovinin and Kiyaeva 2016; Tokovinin 
2020). The closest binaries are impacted by tides. For MS 
stars with convective envelopes below the Kraft break (Teg 
< 6,200 K; Mı < 1.2 Mo), binaries are tidally circular- 
ized below P < 8 days; above the Kraft break, only bina- 
ries within P < 2 days are fully circularized (Meibom and 
Mathieu 2005; Abt 2006; Moe and Di Stefano 2017; Geller 
et al. 2021b; Zanazzi and Wu 2021; Torres et al. 2021). 
Triple Statistics: About 1496 of solar-type primaries 
are in triples and higher order multiples (Tokovinin 2014b, 
see Fig. 1). The closest binaries within P < 3 days are 
found almost exclusively (9696) in triples, whereas only a 
third of binaries with P ~ 20 days have tertiary compan- 
ions (Tokovinin et al. 2006; Laos et al. 2020). Beyond these 
periods, Tokovinin (2014b) concluded that the joint distri- 
bution f(Pin, Pout) of inner periods P, and outer periods 
Pout is consistent with both companions being indepen- 
dently drawn from the same overall lognormal period distri- 
bution for solar-type binaries, but with the added constraint 
of dynamical stability. Various theoretical studies have pro- 
posed slightly different criteria for dynamical stability in 
triples (Mardling and Aarseth 2001; Georgakarakos 2008). 
Based on their sample of solar-type triples with reliable or- 
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TABLE 1 
MULTIPLICITY STATISTICS OF BROWN DWARFS AND MAIN-SEQUENCE STARS 


CBF (%) WBF (%) Ytrunc Ytrunc Ytrunc 
Survey Mı (Mg) N MF (%) THF (%) CF a<lau  a=107-10tau apga ju) a<la a= 1-107 au — a > 107 a 
Fontanive et al. (2018) 0.019- 0.058 47 86 <2 0.08 £ 0.06 8X6 <2 29411 - - - 
Burgasser (2007) 0.05-0.08 162 15443 060.3 0164-004 - - - 
Close et al. (2003) 0.080- 0.095 39 1943 7 «3 0.19 + 0.07% 16 +6 <3 - 
Allen et al. (2007) 0.06-0.15 361 204 «1 0.20 + 0.04 1443 «04 - 
Winters et al. (2019): 20 pc 0.075-0.15 185 1923 22-11 021-3003 16433 05405 31411 - 
Winters et al. (2019): 20 pc 0.15-0.30 336 2312 36410 0.27+0.03 14+2 48412 6+2 - 
Winters et al. (2019): 20 pe 0.3-0.6 350 3042 63-14 0.38 +0.03 1532 1232 1443 - 
EI-Badry et al. (2019): 200 pe 0.1-0.4 5,220 - - - - 
i 0.4-0.6 12,717 - - 
0.6-0.8 9,542 - - 
0.8-1.2 11,588 - - 
ü 12-25 3,238 - - - - - - 
Raghavan et al. (2010): 25 pc 0.75 - 1.25 454 46 +3 1242 0.60 + 0.04 20 +2 2242 49 +6 294 1520.5 02--04 
i 0.75- 1.00 323 4243 - - - - - - 
ie 1.00-1.25 131 5044 - - - - - - 
Tokovinin (2014b): 67 pc 0.85- 1.5 4847 413 1442 0.62 + 0.04 2432 1932 3145 2204 170.5 04404 —0 204 
De Rosa et al. (2014): 75 pc 16-24 435 - - 28 4 - -13404  —22+04 
Murphy et al. (2018) 14-23 2224 - - - - - -11403 
Moe and Kratter (2021) 16-24 - 68 +7 255 0.99 +0.13 3745 - 3248 1333 - - 
Moe and Di Stefano (2017)P 3-5 81-6 3648 128 +0.17 467 40+6 2837 82 05406 —1.0+0.5 
Moe and Di Stefano (2017)P 5-8 89-5 4541 1.55 + 0.24 5448 4937 25:7 63415 03405 — —17405 
Moe and Di Stefano (2017)P 8-17 9344 5715 182-03 61+ 10 55-8 2347 42411 01404 —16+0.5 
Sana et al. 17-50 +4 68+18 21403 70 +11 629 1926 17405 —01£+05 —14+04 


Nore.—All statistics are computed after correcting for incompleteness and removing sys 
we add 4% to their reported resolved binary fractions to account for unresolved spectro 


tems with WD companions. (a): To measure MF and CF from the BD imaging surveys, 
scopic binaries. (b): A compilation of B-type multiplicity surveys, including Abt et al. 


(1990), Shatsky and Tokovinin (2002), and Rizzuto et al. (2013). (c): From the combined Sana et al. (2012) and Sana et al. (2014) O star surveys. 


bital solutions, Tokovinin (2004) measured an empirical re- 
lation: 


5 

(1 in €out)? , 
where €out is the eccentricity of the outer tertiary. 

Tokovinin (2017) examined the relative orbital orienta- 
tions of visually resolved triples. For outer tertiaries beyond 
Gout > 1,000 au, they found an equal number of prograde 
and retrograde orbits with respect to the inner binaries, sug- 
gesting isotropic orientations. Meanwhile, the majority of 
compact solar-type triples within aou; < 50 au have pro- 
grade configurations, demonstrating a high degree of align- 
ment between the inner and outer orbits. Tokovinin (2017) 
also found tentative evidence that triples become more mis- 
aligned with increasing primary mass. Finally, Borkovits 
et al. (2016) identified a large population of very compact 
triples with aout < 10 au. All have mutual inclinations 
within [mutual < 60°, about half of which are nearly copla- 
nar with Jinutual < 20°. We discuss the special implications 
of triple star statistics for formation models in §3.3.3. 


Pout / Pin > (4) 


2.2.2. M-dwarfs 


Frequencies: Building on the historical work of Fischer 
and Marcy (1992), many subsequent studies have employed 
a variety of detection techniques to fill in the multiplicity 
parameter space of the most common stars in the galaxy 
(Close et al. 2003; Gizis et al. 2003; Bouy et al. 2003; 
Daemgen et al. 2007; Law et al. 2008; Bergfors et al. 2010; 
Janson et al. 2012; Dieterich et al. 2012; Janson et al. 2014; 
Ward-Duong et al. 2015). By combining previous samples 
with their own imaging survey, Winters et al. (2019) pro- 
vided the most complete census of the multiplicity statistics 
for a 25-pc volume-limited sample of M-dwarfs primaries. 
We split their sample into three mass intervals and focus 
on the 20-pc subset that is relatively complete (see their 


Figs. 19-20 and results in our Table 1). The multiplicity 
fraction, triple fraction, and median companion separation 
all increase with primary mass within the M-dwarf regime 
(see Figs. 1-3). As shown in Fig. 4, the wide binary frac- 
tion dramatically increases from 1% to 12% across primary 
masses 0.1 - 0.5 Mo. In Fig. 2, we plot the separation dis- 
tribution of early-M (M; = 0.4 -0.6 M) and late-M (Mj = 
0.075 - 0.15 Mo) binaries from the 20-pc subset of Winters 
et al. (2019) along with their best-fit lognormal distributions 
(see Table 2 for parameters). 
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Fig. 3.—: Median separations of all companions (thick black) 
and of inner binaries only (thin green) as a function of pri- 
mary mass. Metal-poor solar-type binaries are skewed toward 
smaller separations (blue) while metal-rich solar-type binaries fa- 
vor longer periods (red). A significant fraction of massive stars 
are in triples, so their inner binary distributions are substantially 


skewed toward shorter separations compared to all companions. 


Winters et al. (2019) excluded both BD and WD com- 
panions from their reported multiplicity statistics but kept 


Offner, Moe, Kratter, Sadavoy, Jensen & Tobin 


Table 2:: LOGNORMAL FIT PARAMETERS 
Mı (Mo) u (au) loga CF 
0.075 - 0.15 4 0.7 0.21 
0.3 - 0.6 25 1.3 0.38 


0.75 - 1.25 40 1.5 0.60 


NOTE - Mean p and dispersion Giog a of lognormal separation dis- 
tributions fiog a scaled to the companion frequency CF integrated 
across —2.0 « log a (au) « 4.5 as defined in Eqn. 3. 


records of such systems in their sample. The fraction of M- 
dwarfs within 20 pc that have observed BD companions is 
only 2%, which increases slightly to 4% after accounting 
for incompleteness. Whether or not BD companions are in- 
cluded, the multiplicity fraction of late-M dwarfs is lower 
than that of early-M dwarfs, which in turn is less than that 
of solar-type primaries. 

Mass ratios: Across intermediate separations, early-M 
binaries follow a roughly uniform mass-ratio distribution, 
similar to that of solar-type binaries. Meanwhile, mid-M 
and late-M binaries are skewed substantially toward equal 
masses, which is not a selection effect due to incomplete- 
ness or exclusion of low-mass companions (Close et al. 
2003; Law et al. 2008; Dieterich et al. 2012; Janson et al. 
2014). Wide companions (a > 100 au) to M-dwarfs are 
skewed toward slightly smaller masses (Ward-Duong et al. 
2015; El-Badry et al. 2019), similar to the separation trend 
Observed in solar-type binaries. 


2.2.8. Brown dwarfs 


Frequencies: Although completeness corrections be- 
come more challenging with decreasing primary mass, sur- 
veys confirm that the MF continues to decline through the 
BD regime (see Table 1 and Fig. 1). The period-integrated 
MF is roughly + 20% for L/early-T dwarfs and only ~ 896 
for late-T/early-Y dwarfs (Reid et al. 2001; Close et al. 
2002; Burgasser et al. 2003; Bouy et al. 2003; Close et al. 
2003; Gizis et al. 2003; Allen et al. 2007; Burgasser 2007; 
Aberasturi et al. 2014; Basri and Reiners 2006; Joergens 
2008; Blake et al. 2010; Hsu et al. 2021; Fontanive et al. 
2018). The vast majority of BD binaries have a — 1 - 20 au. 
The binary fraction within a < 1au is only ~ 496, and only 
a handful of systems have been identified at wide separa- 
tions (Burgasser 2007; Allen et al. 2007; Joergens 2008; 
Radigan et al. 2009; Faherty et al. 2020). The lognormal 
separation distribution of BD binaries narrowly peaks near 
Gay = 3 au, similar to that of late-M binaries (see Fig. 3). 

Mass Ratios: BD binaries are extremely skewed toward 
equal masses, 1.e., Ytrunc = 2 - 3 for L/early-T binaries (Bur- 
gasser et al. 2003; Close et al. 2003; Allen et al. 2007) 
and possibly up to ytrunc © 5 for late-T/early-Y binaries 
(Fontanive et al. 2018). These measurements are not bi- 
ased by incompleteness to smaller mass ratios, nor do they 
exclude objects below the deuterium burning limit. 
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Fig. 4.—: The close (top) and wide (bottom) binary fractions of 
MS stars (thick black) both increase with primary mass but with 
markedly different trends. Metal-poor (orange) solar-type stars 
have a larger CBF than their metal-rich (cyan) counterparts (Moe 
et al. 2019), while the CBF of OB stars is constant with metallicity 
(Moe and Di Stefano 2013). The WBF of solar-type stars is also 
metallicity invariant (Moe et al. 2019; El-Badry and Rix 2019). 
The CBF of Class II/III T Tauri stars (magenta) increases with 
mass like the field MS population (Kounkel et al. 2019). We high- 
light discrepancies in the WBF between the dense Orion Nebula 
Cluster (red; Duchéne et al. 2018) and sparse Taurus star-forming 
region (blue; Kraus et al. 2011), which exhibit a deficit and ex- 
cess, respectively, compared to the field MS. The WBF of Her- 
big Ae/Be stars (green) is consistent with their field MS analogs 
(Kouwenhoven et al. 2005; Baines et al. 2006), but the CBF of 
Herbig Ae/Be stars appears lower (Corporon and Lagrange 1999; 
Apai et al. 2007; Sana et al. 2017); “?” indicates this could arise 
from a selection bias whereby close binaries shorten disk lifetimes. 


2.2.4. A stars 


Frequencies: Moving toward earlier spectral types, we 
find that the trend of an increasing MF with primary mass 
continues. The wide binary fraction of A-dwarfs is slightly 
larger than the solar-type value (De Rosa et al. 2014). 
Early spectroscopic binary studies of A stars focused on 
their magnetic and chemical peculiar properties (Abt 1965, 
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1983). There has been no recent systematic A-dwarf spec- 
troscopic or eclipsing binary survey, but the upcoming Gaia 
full third data release will soon remedy this ($6.2). Instead, 
close binary properties of A stars have been inferred from 
examining Kepler phase modulations of ó Scuti stars, which 
are pulsating late-A/early-F MS stars (Murphy et al. 2018). 
They found an A-dwarf binary fraction of 14% across 0.6 - 
3.5 au, double the G-dwarf value across the same separation 
interval. Since the data remain incomplete, we use an in- 
terpolation from Moe and Kratter (2021) to derive the full 
bias-corrected multiplicity statistics for A-dwarfs across the 
entire separation range (see Table 1 and Fig. 1). 

Mass ratios: Across intermediate separations, A-dwarf 
binaries peak near g = 0.3 and exhibit a very small but statis- 
tically significant excess twin fraction (De Rosa et al. 2014; 
Moe and Di Stefano 2017; Murphy et al. 2018). Meanwhile, 
wide A-type binaries are skewed toward even smaller mass 
ratios but are still slightly top-heavy compared to random 
pairings drawn from the IMF (De Rosa et al. 2014; Moe 
and Di Stefano 2017). 


2.2.5. B stars 


Frequencies: Building on the groundbreaking work by 
Abt et al. (1990), surveys of B-stars show binary fractions 
that continue to rise with primary mass across all orbital 
separations (Duchéne et al. 2001; Shatsky and Tokovinin 
2002; Roberts et al. 2007; Rizzuto et al. 2013; Moe and Di 
Stefano 2015b; Caballero-Nieves et al. 2020; Banyard et al. 
2022). Moe and Di Stefano (2017) combined the data from 
spectroscopy, eclipsing binaries, long-baseline interferom- 
etry, and adaptive optics to compute the statistics reported 
in Table 1. Both the masses and companion frequencies of 
B stars span a large interval, increasing from CF = 1.1 for 
3 Mo primaries to 1.9 for 16 Mọ primaries (see Fig. 1). 

The overall separation distribution of companions to 
mid-B stars roughly follows Opik’s law, i.e., uniform in 
log a. (Öpik 1924; Kouwenhoven et al. 2007; Kobulnicky 
and Fryer 2007; Moe and Di Stefano 2017). Fig. 2 displays 
fioga for q > 0.1 companions to 8 Me mid-B primaries as 
reported in Moe and Di Stefano (2017). Early-B inner bina- 
ries favor shorter separations with a median of à; 7 6 au 
(Rizzuto et al. 2013; Moe and Di Stefano 2013, see Fig. 3). 

Mass ratios: B-type binaries within a < 1 au follow a 
uniform mass-ratio distribution (Abt et al. 1990; Kobulnicky 
et al. 2014). Companions to B-type stars across intermedi- 
ate separations of a = 1- 100 au peak near q ~ 0.3 (Riz- 
zuto et al. 2013; Gullikson et al. 2016; Moe and Di Stefano 
2017), similar to their A-type counterparts. Wide compan- 
ions to B-type stars are substantially skewed toward small 
mass ratios (7trunc * —2) but with a slight flattening below 
q S; 0.3 (Abt et al. 1990; Shatsky and Tokovinin 2002; Moe 
and Di Stefano 2017). 


2.2.6. O stars 


Frequencies: The total multiplicity fraction of O stars 
in clusters is MF >90% and the companion frequency is 
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CF = 2.1 +0.3 above q > 0.1 (see Figs. 1-4). Hence, 
the majority of O stars are in triples and higher ordered 
multiples (THF = 68% + 18%). As with solar-type stars, 
these statistics derive from a wide range of survey methods: 
spectroscopy, interferometry, sparse aperture masking, AO, 
lucky imaging, and common proper motion (Turner et al. 
2008; Mason et al. 2009; Sana et al. 2012; Chini et al. 
2012; Kobulnicky et al. 2014; Sota et al. 2014; Sana et al. 
2014; Aldoretta et al. 2015; Moe and Di Stefano 2017; Maíz 
Apellániz et al. 2019; Caballero-Nieves et al. 2020). Given 
the high fraction of triples, we display both the total and 
inner-binary separation distributions in Fig. 2. 

The statistics above describe O stars in young clusters. 
In contrast, the ~ 20% of O stars that are runaway and/or in 
the field exhibit a lower binary fraction (Mason et al. 2009; 
Chini et al. 2012; Lamb et al. 2016). The multiplicity statis- 
tics of older field and runaway O stars have been altered 
by dynamical effects and/or binary evolution (Hoogerwerf 
et al. 2001). Hence, we consider only cluster O stars, which 
have multiplicity properties more similar to their primordial 
distributions. 

Mass Ratios: Spectroscopic and eclipsing binary sur- 
veys of O stars reveal a fairly uniform mass-ratio distribu- 
tion within a « 0.5 au but with a small excess twin fraction 
(Pinsonneault and Stanek 2006; Sana et al. 2012; Kobul- 
nicky et al. 2014; Moe and Di Stefano 2017; Shenar et al. 
2022). Nonetheless, the excess twin fraction among O-type 
binaries within a « 0.5 au is only 10%, which is lower than 
the 3096 excess twin fraction measured among solar-type 
binaries across the same separation interval (Moe and Di 
Stefano 2017). There is no excess twin fraction among O- 
type binaries beyond a > 0.5 au (Sana et al. 2012). Similar 
to early-B binaries, the mass-ratio distribution of O-type bi- 
naries beyond a > 100 au is skewed steeply toward small 
mass ratios but with a flattening below q « 0.3 (Sana et al. 
2014; Moe and Di Stefano 2017). 


2.2.7. Summary of observed MS trends 


The multiplicity fraction increases monotonically with 
primary mass from MF ~ 20% for BDs and late-M dwarfs 
to 225096 for solar-type stars to MF > 90% for OB stars 
(see Fig. 1). The triple fraction increases even more dra- 
matically from THF ~ 2% for late-M dwarfs to 14% for 
FGK dwarfs to nearly 7096 for O stars. Similarly, the com- 
panion frequency is CF = 0.2 for low-mass stars while OB 
stars have CF ~ 2.0 companions per primary on average. 

Figs. 2-4 illustrate the nuances of the separation distri- 
bution with respect to primary mass. BD, M-dwarf, and 
FGK-dwarf binaries follow lognormal separation distribu- 
tions, where the peak increases from Gay = 3 au for BD and 
late-M dwarfs to à4j = 40 au for solar-type stars. Mean- 
while, companions to OBA stars roughly follow Opik's law, 
albeit inner binary companions to early-B and especially 
O stars are skewed toward shorter separations, e.g., di, © 
2 au. Interestingly, the close binary fraction is nearly con- 
stant at CBF = 15% across M, = 0.05-0.8 Mo but then 
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rapidly increases with primary mass, reaching 7096 for O 
stars. Meanwhile, the wide binary fraction increases dra- 
matically from WBF = 1% for late-M stars to 12% for early- 
M stars, and then continues to moderately increase up to 
60% for OB stars. The different trends in the close versus 
wide binary fractions can help constrain the various binary 
formation processes that operate on different scales. 

The binary mass-ratio distribution becomes skewed 
toward more unequal masses with increasing primary 
mass and increasing orbital separation. It is firmly es- 
tablished that BD and late-M binaries intrinsically favor 
equal masses. We note that some binary populations fa- 
vor mass-ratio distributions that deviate from a power-law 
model: close and intermediate-period solar-type binaries 
follow uniform distributions with a 3096 and 1096 excess 
twin fraction, respectively (see $2.2.1), very close O-type 
binaries follow a uniform distribution with a 10% excess 
twin fraction (82.2.6), and intermediate-period compan- 
ions to BA primaries peak near q = 0.3 ($2.2.4-2.2.5). Fi- 
nally, wide companions to OBA primaries match the slope 
Ytrunc = —2.35 expected from random pairings drawn from 
a Salpeter IMF, but only across the truncated interval q = 
0.4 - 1.0. The mass-ratio distribution flattens below q « 0.3, 
so the overall distribution is slightly top-heavy compared to 
random pairings from the IMF (see $2.2.4-2.2.6). Observed 
mass-ratio distributions demonstrate that the component 
masses in close and intermediate-period binaries coevolved 
during their prior fragmentation and accretion phases, while 
the components in wide binaries formed mostly, but not 
fully, independently (see $3). 


2.3. Metallicity Dependence 


Since PPVI, the field has seen a substantial change in 
our understanding of the metallicity dependence of binary 
statistics. Early surveys of metal-poor halo stars measured 
a spectroscopic binary fraction that is consistent with the 
disk population (Latham et al. 2002; Carney et al. 2005). 
However, it is more difficult to detect spectroscopic binary 
companions to metal-poor stars, which have weaker absorp- 
tion lines, and so the bias-corrected close binary fraction of 
halo stars is definitively larger (Moe et al. 2019). There are 
now several surveys indicating that the close binary fraction 
of solar-type stars decreases significantly with metallicity 
(Grether and Lineweaver 2007; Rastegaev 2010; Raghavan 
et al. 2010; Badenes et al. 2018; Moe et al. 2019; Mazzola 
et al. 2020). For example, both APOGEE spectroscopic bi- 
naries and Kepler eclipsing binaries reveal a close binary 
fraction that decreases by a factor of four across —1.0 « 
[Fe/H] « 0.5 (Moe et al. 2019, see Fig. 4). Mazzola et al. 
(2020) measured an even steeper anti-correlation between 
the close binary fraction and chemical abundances of o el- 
ements like C, O, Mg, and Si. 

Meanwhile, the wide binary fraction across a — 200- 
1,000 au appears to be metallicity invariant (Moe et al. 
2019; El-Badry and Rix 2019). El-Badry and Rix (2019) 
utilized their catalog of common-proper-motion binaries 
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from Gaia to demonstrate that a strong metallicity depen- 
dence emerges only within a < 200 au. They measured a 
factor of three decrease in the fraction of solar-type binaries 
near a = 50 au across —1.0 < [Fe/H] < 0.5, nearly the fac- 
tor of four variation observed within a < 10 au across the 
same metallicity interval. Figs. 2 & 3 show that metal-poor 
solar-type binaries skew toward shorter separations while 
metal-rich binaries favor long periods. These observations 
indicate metallicity plays a crucial role in binary fragmen- 
tation at high densities and small separations (see §3.4.5). 

Finally, the small population of very wide solar-type 
binaries across a = 1,000-10,000 au exhibits a non- 
monotonic metallicity dependence in their companion fre- 
quency, peaking near [Fe/H] = 0.0 (Hwang et al. 2021). 
The deficit of extremely wide companions (a > 5,000 au) 
to older, metal-poor halo stars is at least partially due to 
dynamical disruption via galactic tides and gravitational 
perturbations. However, Hwang et al. (2021) argued that 
disruptions alone cannot explain the metallicity dependence 
observed in binaries across a = 1,000 - 5,000 au, which have 
larger binding energies. They instead concluded that metal- 
poor halo stars initially formed in denser clusters, result- 
ing in a smaller very wide binary fraction (see also §2.4). 
Meanwhile, Hwang et al. (2021) offered multiple theories 
for the smaller very wide binary fraction of [Fe/H] = 0.5 
disk stars, including radial migration within the Milky Way 
and/or dynamical unfolding of initially unstable triples. 

Unlike solar-type stars, the close binary fraction of mas- 
sive stars does not decrease with metallicity. Both eclipsing 
and spectroscopic binary surveys of OB stars in the metal- 
poor Magellanic Clouds reveal a close binary fraction that 
matches their Milky Way counterparts at solar-metallicity 
(Moe and Di Stefano 2013; Sana et al. 2013; Dunstall et al. 
2015, see Fig. 4). Moreover, the period and mass-ratio dis- 
tributions of close OB binaries in the Magellanic Clouds are 
consistent with the Milky Way distributions (Moe and Di 
Stefano 2013; Almeida et al. 2017; Villasenor et al. 2021). 
Recent surveys of massive metal-poor stars in older envi- 
ronments indicate a slightly lower binary fraction (Boden- 
steiner et al. 2021; Neugent 2021), but this may be due to 
binary evolution and/or other selection effects. 


2.4. Pre-MS Phase 


2.4.1. T Tauri stars 


Spectroscopic surveys of Class II/III T Tauri stars reveal 
a close binary fraction that is consistent with their field MS 
counterparts (Mathieu 1994; Melo 2003; Elliott et al. 2014; 
Prato 2007; Kounkel et al. 2019). Specifically, Kounkel 
et al. (2019) measured a bias-corrected close binary frac- 
tion that increases from CBF = 12% for mid-M T Tauri stars 
to CBF = 33% for early-F T Tauri stars (see Fig. 4). They 
also found that close T Tauri binaries follow the same short- 
period tail of the lognormal separation distribution as field 
MS binaries across a = 0.1- 10 au (see Fig. 2). These ob- 
servations demonstrate that the primary mass dependence 
and separation distribution of close binaries are largely set 
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by the beginning of the T Tauri stage. 

Two recent spectroscopic surveys suggest that the close 
binary fraction decreases with age across T = 2- 100 Myr 
(Jaehnig et al. 2017; Zuifiiga-Ferndndez et al. 2021), 
but we argue these results are due to selection biases. 
Based on APOGEE observations, Jaehnig et al. (2017) 
reported completeness-corrected CBF = 1296-2096 for 
M, = 0.5 Mo pre-MS stars in five different star-forming 
groups, which is consistent with the early-M dwarf field 
value. In contrast, they measured only CBF ~ 4% for 
low-mass MS stars in the ~ 100 Myr old Pleiades clus- 
ter. This discrepancy likely arises because the low-mass 
stars in the Pleiades observed by APOGEE are concen- 
trated toward the periphery of the cluster, which has a lower 
binary fraction due to mass segregation (Raboud and Mer- 
milliod 1998). Moreover, the close binary fraction of FGK 
MS stars in the Pleiades measured from more robust long- 
term radial velocity monitoring is CBF = 25% + 396 (Tor- 
res et al. 2021), similar to the corresponding field value 
(see 82.2.1). Züfiiga-Fernández et al. (2021) recently mea- 
sured an even larger solar-type close binary fraction of CBF 
= 30% + 6% averaged across three young moving groups 
with 7 « 20 Myr, including Beta Pic, and measured only 
CBF z 1096 - 1596 in slightly older associations with 7 — 
20-100 Myr. However, of their ten spectroscopic binaries 
in Beta Pic, six have slow rotation velocities and/or peri- 
ods compared to their single star analogs, suggesting they 
are older non-members (only two have short orbital peri- 
ods where tides could have reduced the spin). Moreover, 
four of the ten have > 20% probabilities of not being Beta 
Pic members, while only two of the remaining 32 members 
have such low reliabilities. After removing likely non- 
member spectroscopic binaries, the close binary fraction in 
young comoving groups matches the field value. 

Meanwhile, AO and speckle imaging surveys of T Tauri 
stars in star-forming environments like Taurus, Chamaeleon, 
Ophiuchus, and Scorpius reveal an excess of companions 
across a = 10 - 200 au (Ghez et al. 1993; Reipurth and Zin- 
necker 1993; Leinert et al. 1993; Ghez et al. 1997; Kraus 
et al. 2011; Tokovinin and Bricerio 2020), which is discon- 
tinuous with the spectroscopic binary results (see Fig. 2). It 
was originally assumed that only low-density environments 
exhibited such an excess while high-density environments 
would have a deficit, and therefore the field population 
derived from a combination of these environments. How- 
ever, Duchéne et al. (2018) found tentative evidence for 
the same excess of companions across a = 10-60 au in 
the dense Orion Nebula Cluster (ONC). Notably De Furio 
et al. (2019) found no excess in this separation range for 
low-mass M-dwarfs in Orion. Thus the discontinuity and 
excess near a — 10 au in most star-forming regions poses a 
mystery and warrants further study. 

Across wide separations of a = 100-10,000 au, the 
binary fraction is sensitive to the stellar density of the 
environment. Fig. 2 shows that the dense Orion Neb- 
ula Cluster (ONC) exhibits a deficit of wide companions 
(Scally et al. 1999; Kóhler et al. 2006; Reipurth et al. 
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2007; Duchéne et al. 2018; Jerabkova et al. 2019). En- 
vironments with intermediate density such as Upper Scor- 
pius and the Orion OBI association are consistent with the 
field (Brandner et al. 1996; Kraus et al. 2008; Kounkel 
et al. 2016; Tokovinin and Briceño 2020; Tokovinin et al. 
2020). Finally, low-density star-forming regions like Tau- 
rus and Chamealeon exhibit a slight excess of wide com- 
panions (Ghez et al. 1993; Reipurth and Zinnecker 1993; 
Leinert et al. 1993; Ghez et al. 1997; Kohler and Leinert 
1998; Connelley et al. 2008; Kraus et al. 2011; Joncour 
et al. 2017). These observations suggest that the majority of 
wide binaries are dynamically disrupted in extremely dense 
regions, and that the field MS population derives from a 
mixed composition of low, intermediate, and high-density 
environments. Specifically, for every solar-type primary 
born in a dense ONC environment, a counterpart must be 
born in a sparse Taurus-like association such that the aver- 
age of the star-forming environments yields a wide binary 
fraction similar to Upper Scorpius and the field MS popula- 
tion. 

However, within the same star-forming region, there 
is tentative evidence for the opposite trend between the 
wide binary fraction and the immediately surrounding stel- 
lar density. For example, Kounkel et al. (2016) imaged 129 
Class I protostars and 197 Class II/III T Tauri stars across 
the Orion Molecular Cloud in the near-IR with HST and 
the InfraRed Telescope Facility. Across a = 100 - 1,000 au, 
they found that protostars and T Tauri stars have very simi- 
lar wide binary fractions, which are consistent with the field 
MS population. Surprisingly, Kounkel et al. (2016) found 
that young stars in sub-regions of higher spatial stellar den- 
sity exhibited a slightly larger wide binary fraction. Recent 
ALMA observations of Class 0/I protostars in Orion seem 
to confirm this result (Tobin et al. 2022; see 84.2.1). 

There are noticeable differences between close versus 
wide T Tauri binaries. Similar to their MS counterparts, 
T Tauri binaries within a « 200 au follow a uniform 
mass-ratio distribution with a small excess twin fraction, 
while wider companions are skewed toward smaller masses 
(Kohler and Leinert 1998; Tokovinin and Bricefio 2020). As 
discussed in §5.1.1, binaries within a < 50 au truncate the 
disk masses, radii, fluxes, and lifetimes compared to single 
stars and wider binaries (Jensen et al. 1996; Harris et al. 
2012; Kraus et al. 2012; Cheetham et al. 2015). Hence, a 
sample of T Tauri stars selected according to their large disk 
masses will be systematically biased against close binaries. 
For example, Kounkel et al. (2019) measured the fraction of 
Class II and III T Tauri stars that are in double-lined spec- 
troscopic binaries (SB2s), which have large mass ratios q > 
0.6. Compared to the field MS population, they found that 
Class II stars exhibit a deficit of SB2s while Class III stars 
have a surplus. Kounkel et al. (2019) concluded that close 
T Tauri binaries with large mass ratios quickly consume or 
disrupt their disks, thereby accelerating the Class II phase 
and extending the duration of the Class III phase compared 
to single stars and wider binaries. 
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2.4.2. Herbig Ae/Be stars and Massive YSOs 


According to various high-resolution imaging surveys, 
the wide binary fraction of Herbig Ae/Be stars and mas- 
sive YSOs in young star-forming environments is WBF 
= 30% - 60%, increasing with stellar mass (Kouwenhoven 
et al. 2005; Baines et al. 2006; Wheelwright et al. 2010; 
Gravity Collaboration et al. 2018; Pomohaci et al. 2019). 
These values are consistent with their MS analogs in the 
field and slightly older clusters (see Fig. 4). Wide compan- 
ions to OB stars in the Trapezium Cluster of the ONC skew 
toward small mass ratios q ~ 0.2, similar to their counter- 
parts in slightly older clusters (Gravity Collaboration et al. 
2018). The samples of wide Herbig Ae/Be binaries are 
not sufficiently large to detect a dependence on surround- 
ing stellar density as measured for T Tauri stars. 

Meanwhile, the close binary fraction of Herbig Ae/Be 
stars and massive YSOs inferred from spectroscopy is only 
CBF = 10% -30% (Corporon and Lagrange 1999; Apai 
et al. 2007; Sana et al. 2017), considerably lower than 
the corresponding MS fraction. One interpretation is that 
1 Myr old massive stars do not have any close companions, 
and that the companions at intermediate separations harden 
on Myr timescales (see also Ramírez-Tannus et al. 2021). 
An equally compelling explanation is that close compan- 
ions to massive pre-MS stars shorten disk lifetimes, similar 
to their T Tauri analogs, and therefore the subset of mas- 
sive stars within a coeval cluster that retains a disk is less 
likely to have close companions. Resolving this ambiguity 
requires a systematic survey for close binary companions 
to both massive YSOs with disks and normal OBA stars 
without disks within the same cluster (Bordier et al. 2022). 
Larger samples of spectroscopic and eclipsing OB binaries 
in young star-forming environments within the Large Mag- 
ellanic Cloud reveal a sizeable fraction with close compan- 
ions (Sana et al. 2013; Moe and Di Stefano 20152). For ex- 
ample, Moe and Di Stefano (20152) discovered 18 early-B 
MS stars with 1-2 Mọ pre-MS eclipsing binary compan- 
ions across P = 3-8 days, a few of which are as young as 
0.6 - 1.0 Myr according to the measured temperatures and 
inflated radii of the pre-MS companions still on the Hyashi 
track. 


2.5. Overall Multiplicity Statistics 


Given the substantial excess of wide companions to 
YSOs compared to field stars (see $2.4 and $4), most star 
systems are born as multiples. However, subsequent dy- 
namical disruptions reduce the overall multiplicity fraction, 
especially for low-mass binaries that have lower binding en- 
ergies (Kroupa 1995). As emphasized in Lada (2006), most 
MS field stars are M-dwarfs and most M-dwarf field stars 
are single, and thus most MS star systems are single. By 
weighting our multiplicity statistics with the IMF of pri- 
mary stars (Kroupa et al. 2013), we compute that only 3596 
of zero-age MS star systems are multiple, consistent with 
the conclusion of Lada (2006). Nonetheless, most MS stars 
are members of multiples after considering the fact that 
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multiples contain two or more stars. For example, given 
100 M-dwarfs, only 49 are single, 16 are paired with 16 
other M-dwarfs, 12 are binary companions to more mas- 
sive BAFGK primaries, and 7 are tertiaries in triples and 
higher-ordered multiples. By weighting our statistics with 
the IMF, we conclude that 58% of MS stars are members of 
multiples. 


3. MODELS FOR MULTIPLE STAR FORMATION 


A successful theory for multiple formation must both 
reproduce observed multiplicity statistics ($2) and model 
young, individual sources ($4 and $5). To date, there is no 
comprehensive, predictive theory for multiple star forma- 
tion. Models have instead tackled the origin of multiplicity 
at different spatial scales and explored the influence of dif- 
ferent physical processes. 

By definition, multiple formation requires two or more 
events of gravitational collapse that ultimately form a 
bound stellar system. We define the process by which 
self-gravitating objects develop substructures that evolve 
and collapse independently as fragmentation. Hoyle (1953) 
first proposed star formation through the hierarchical frag- 
mentation of collapsing gas. However, this “classic” star 
formation paradigm, described only by the interplay of 
gravity and thermal pressure, produces smooth, isotropic 
collapse (e.g., Shu 1977). Some of the first numerical star- 
formation calculations demonstrated that multiple fragmen- 
tation events are not possible under these conditions alone 
(Larson 1972). 

Observations and theoretical studies have reinforced that 
this simple picture is a poor representation of how most 
stars form. Instead, star formation is regulated not only 
by gravity and thermal pressure but also turbulence and 
magnetic fields, which introduce rotation and asymmetry, 
crucial ingredients that facilitate additional collapse events 
within dense cores, filaments and accretion disks. The high 
degree of clustering in many star-forming regions also sug- 
gests that interactions between stars are important for the 
evolution of multiple systems. 

Current ideas for multiple formation can be divided into 
three main categories: theories in which multiples form via 
fragmentation of a core or filament ($3.1), via fragmenta- 
tion of a massive accretion disk ($3.2) or through dynamical 
interactions (§3.3). This third mode can also rearrange the 
hierarchy and multiplicity of systems formed via the prior 
fragmentation channels. Figure 5 summarizes these mecha- 
nisms and provides examples from observations and numer- 
ical simulations. $3.4 reviews the landscape of state-of-the- 
art calculations, approaches the origin of stellar multiplicity 
holistically and discusses our current understanding of how 
physics and environment shape multiplicity statistics. 


3.1. 


Stars form in dense cores, most of which are observed 
to be embedded within ~ 0.1 pc filaments in nearby star- 
forming regions (Könyves et al. 2015). High-resolution ob- 


Core and Filament Fragmentation 
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Fig. 5.—: Summary of mechanisms for multiple formation. Top: Model and approximate range of time and length scales for each 
process. Middle: Proposed observational examples. From left to right: B5 in Perseus (Pineda et al. 2015), SMIN in Ophiuchus (Kirk 
et al. 2017), L1448 IRS3B in Perseus (Reynolds et al. 2021) and RW Aur (Rodriguez et al. 2018). Bottom: Examples from numerical 
simulations. From left to right: Guszejnov et al. (2021) Offner et al. (2016), Bate (2018), and Muñoz et al. (2015). 


servations have discovered even smaller filaments within 
cores that host forming stars (Pineda et al. 2015, see 
Fig. 5). Core and filament fragmentation models posit 
that stellar multiples arise from over-densities that develop 
and collapse within these parent structures, producing ini- 
tially widely separated systems (= 500 au). Note that the 
terms core and filament represent morphological descrip- 
tions rather than physical definitions, and here, we define 
over-densities with small aspect ratios (r;/rp < 3/1) as 
cores and structures with higher aspect ratios as filaments. 
However, we stress that these terms likely represent two 
limits on a continuous spectrum of initial gas morphologies 
and physical conditions rather than wholly distinct classes 
of objects. Below we discuss different processes that regu- 
late the fragmentation of these structures and the predicted 
signatures of their formation. 


3.1.1. 


Two proposed drivers of core fragmentation are rotation 
and turbulence, which produce density and velocity asym- 
metries. Velocity gradients observed in early observations 
of dense cores (e.g., Goodman et al. 1993) highlighted the 
importance of angular momentum in core evolution. Early 
numerical simulations including solid body rotation demon- 
strated that rapidly rotating, collapsing cores are prone to 
fragmentation, thereby leading to binary formation (Larson 
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1972). Inutsuka and Miyama (1992) quantified the frag- 
mentation criterion for rotating cores in terms of the initial 
ratio of thermal to gravitational energy, o, and ratio of ro- 
tational to gravitational energy, Brot, where ovi; fioc > 0.12 
produces fragmentation and 

) Gf)! 


Fragmentation also requires a<0.5, since thermal pres- 
sure may otherwise prevent warm cores from fragmenting 
(Tsuribe and Inutsuka 1999). These criteria are straight- 
forward, but they fail to encapsulate the complexity of real 
cores and thus are not readily extended to construct a gen- 
eral theory of binary formation. 

While recent high-resolution core observations validate 
the prevalence of velocity gradients, current interpretation 
of their physical meaning is nuanced. Gradients may in- 
stead signify core formation via converging flows (Chen 
et al. 2020) or represent the largest turbulent fluctuation 
in the core (Goodwin et al. 2004). The latter explanation 
suggests that turbulence, either directly or indirectly, rather 
than rotation produces multiple star formation (see $3.1.3 
for more discussion). 

Most current star-formation theories appeal to turbu- 
lence to trigger the growth of structure, including the forma- 
tion of cores and filaments, within molecular clouds. Sev- 
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eral groups have developed rigorous, semi-analytic frame- 
works for star formation through turbulent fragmentation, 
in which the statistical properties of turbulence, such as the 
power spectrum and density probability distribution func- 
tion, determine core masses and regulate gravitational col- 
lapse (see review by Lee et al. 2020b). Only one model, 
which we discuss here, makes quantitative predictions for 
stellar multiplicity. Guszejnov and Hopkins (2015) use a 
grid of perturbations to hierarchically compute the frag- 
mentation of a molecular cloud assuming that turbulence 
both seeds collapse and provides pressure support. This 
method preserves the spatial relationships of collapsing ob- 
jects, such that density fluctuations within cores naturally 
lead to fragmentation and multiple formation (Guszejnov 
et al. 2017). Sub-regions collapse when the local gravita- 
tional energy exceeds the combined turbulent and thermal 
energy. This may be phrased as a critical density threshold: 


) [eet "| (6) 


where po, Ro and M are the initial density, radius, and 
Mach number of the parent cloud, respectively. This for- 
malism assumes turbulent energy scales as E(R)ocR? and 
the cloud is isothermal. For supersonic turbulence p—2 
while for subsonic turbulence p—5/3, where the transition 
occurs at the sonic scale, R,. For a virialized molecular 
cloud, where R,~0.65pc and R«0.1pc is the fragment 
radius, we can write (Guszejnov and Hopkins 2015): 
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where the mass of a 10?au radius fragment is M [^ 
4/37 Perit R$2:0.07 M for typical cloud parameters. Gusze- 
jnov et al. (2017) used this formalism to directly predict 
multiplicity properties finding relatively good agreement 
with observed statistics despite their simplistic treatment of 
turbulence (see $3.4 for more discussion). 

Filamentary structure is a natural by-product of turbu- 
lence. Dense filaments are prone to sub-fragmentation at 
regular internals like “beads on a string,” which provides a 
mechanism to form widely-spaced multiples if filaments are 
sufficiently dense and compact (see Fig. 5). An isothermal 
filament collapses if its mass per unit length, Mj, exceeds a 
critical value. A filament threaded laterally by a magnetic 
field, Bi , will become unstable when Mj exceeds (Inutsuka 
and Miyama 1997; Tomisaka 2014): 
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where Bmag=87c3p/ B} is the ratio of thermal to magnetic 
pressure. Several authors have developed semi-analytic 
models deriving the spectrum of fragment masses (i.e., 
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the core mass function) from filament properties (Inutsuka 
2001; André et al. 2019). We can apply these arguments to 
derive expected relationships for multiple formation from 
filament fragmentation. 

In order for density perturbations to grow, analytic ar- 
guments suggest that substructures within filaments must 
have separations of at least 2.5 times the filament full- 
width half maximum, Wtwnm (Fischera and Martin 2012). 
This suggests that fragments that are sufficiently closely 
spaced to remain bound must form in filaments with widths 
of <0.1 pc/2.5~9x10% au. Assuming turbulent support is 
negligible, a filament fragment will collapse if its mass ex- 
ceeds the critical Bonnor-Ebert mass (André et al. 2019): 


4 
C 
MBE, tn ™1.3 o2 (9) 
= 
A T M f Myei Wewhm 
ell) (aoi) (39222) Mo, 


where the filament surface density, “gi, is the ratio of 
the filament mass-per-length to the width of the filament 
X&a1— Miine/Wewhm- Equation 9 suggests that the forma- 
tion of closely-spaced fragments with small masses requires 
a high filament mass-per-length. This implies that strong 
magnetic fields facilitate the formation of dense, narrow fil- 
aments and thus play a critical role in multiples formed via 
filament fragmentation. 

While these expressions are relatively simple and pro- 
vide intuition for the origin of multiplicity in different the- 
oretical frameworks, they make a variety of simplifying as- 
sumptions. Consequently, using these to accurately predict 
future multiplicity from the properties of a particular ob- 
served filament or core is not straight-forward. Moreover, 
debate in the star-formation community continues on the 
relative roles of rotation, turbulence and magnetic fields on 
gas structure formation and evolution. In $3.4, we exam- 
ine the impact of non-linear physical effects in star-cluster 
formation calculations and discuss the impact of different 
processes on multiple formation across all scales. 


3.1.2. Predicted signatures 


Multiple systems formed from core and filament frag- 
mentation have various characteristics imprinted during for- 
mation that shape their statistical distributions and may 
point to their origins long after the natal gas is dispersed. 

Gravitational fragmentation of an isothermal, turbulent 
cloud is scale free and thus requires additional physics to set 
the minimum fragmentation scale (Guszejnov et al. 2020). 
Current theoretical models appeal to either angular momen- 
tum, magnetic support or tidal forces to set this minimum 
scale (Guszejnov et al. 2017; Haugbglle et al. 2018; Lee and 
Hennebelle 2018); all of these effectively impose a lower 
limit of ~10? au. Fragments that collapse with separa- 
tions 20.1 pc are unlikely to be initially bound and will 
naturally drift apart. These two constraints limit the ini- 
tial separations of multiples formed from turbulent frag- 
mentation to ~ 10? au - 0.1 pc. This separation range is 
consistent with those of multiplies that form in magneto- 
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hydrodynamic simulations of turbulent clouds and cores 
(Offner et al. 2016; Lee et al. 2019; Kuffmeier et al. 2019). 
Arc- and bridge-like gas structures, produced by formation 
in a turbulent environment, may connect the pairs and facil- 
itate accretion (Kuffmeier et al. 2019). Wide embedded and 
accreting, multiples exist at such separations for a relatively 
short time period, migrating to <10? au within ~ 10° years 
(Offner et al. 2016; Lee et al. 2019, see $3.3). 

Turbulent fragmentation is only efficient when there is a 
significant reservoir of gas, i.e., during the embedded phase. 
Consequently, multiples in such systems, assuming they 
form from the same parent structure, are expected to have 
similar ages (AtZtg 3 x 10? (ng, /10*cm ?) 9 yr). 
Once collapse commences, gravity curtails the growth of 
perturbations, limiting the number of objects that can form 
by this process. Consequently, turbulent fragmentation may 
only produce two or three members within a core (Offner 
et al. 2016; Guszejnov et al. 2017; Lee et al. 2019). Since 
collapse events happen independently and occur while there 
is a significant gas reservoir, the formation of very low-mass 
companions (< 0.08M5) and systems with high mass ratios 
are expected to be rare (Fisher 2004). 

Due to their wide initial separations, multiples formed 
from turbulent fragmentation accrete gas with different net 
angular momentum. This frequently produces misaligned 
stellar spins, accretion disks and protostellar outflows (Bate 
2018; Offner et al. 2016; Lee et al. 2019). Stellar spins 
and outflows of wide-separation protostellar pairs formed in 
magneto-hydrodynamic (MHD) simulations are consistent 
with a random distribution (Offner et al. 2016; Lee et al. 
2019). The misalignment persists to the end of the cal- 
culation, suggesting that misaligned spins are a signature 
of this formation mechanism that may be observable after 
the birth cloud disperses (see $4.3.4). The relative orienta- 
tions of protostellar outflows, which are observable before 
significant dynamical evolution occurs, provide a means to 
directly probe the initial angular momentum direction (see 
$4.3.4). However, outflow angles can be difficult to mea- 
sure for close pairs and in clustered regions; selection ef- 
fects due to the difficulty of distinguishing aligned proto- 
stellar outflows may cause wide pairs to appear preferen- 
tially anti-aligned rather than randomly distributed (Offner 
et al. 2016). 


3.1.3. Rotation, turbulence, and angular momentum in- 
heritance 


Several authors have drawn a direct line between the 
angular momentum generated by turbulence and multiplic- 
ity properties. Fisher (2004) developed a semi-empirical 
model relating randomly drawn turbulent velocity fields to 
binary angular momentum, which in turn set binary pe- 
riods and separations. This approach, which depends on 
several input parameters such as the core-to-star formation 
efficiency correctly predicts that more equal-mass binaries 
have smaller binary periods and that lower-mass systems 
have smaller typical binary separations. Jumper and Fisher 
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(2013) extended this model to explore the origin of the 
"brown dwarf desert," in which solar mass stars have rel- 
atively few brown dwarf companions at short periods (see 
§2.2.1). They showed that low-mass brown dwarfs have 
smaller mean separations, where the small fraction of ob- 
served wide-separation brown dwarf systems are explained 
by turbulent stochasticity. 

Groups carrying out hydrodynamic calculations of tur- 
bulent molecular clouds both with and without magnetic 
fields have shown that turbulence naturally produces ve- 
locity gradients with the magnitude and distribution seen 
in observations (e.g., Chen et al. 2020). However, nu- 
merical studies of fragmentation in collapsing, magnetized 
cores show poor correlation between the initial degree of 
rotation and later fragmentation (Kuznetsova et al. 2020). 
Kuznetsova et al. (2020) instead propose that dynamical 
torques are more directly responsible for core fragmenta- 
tion. Ultimately, multi-physics simulations highlight the 
difficulty of applying simple theoretical models to predict 
the complex and nonlinear outcomes of star formation. 

Non-ideal magnetic effects, which influence angular mo- 
mentum transport on both core and filament scales, further 
confuse the relationship between initial gas properties and 
fragmentation (Zhao et al. 2020). The origin and distribu- 
tion of angular momentum is in turn crucial for the forma- 
tion and properties of protostellar disks, which may also 
host binary formation as we discuss in the next section. In 
§3.4, we present a more detailed discussion of multiplicity 
derived from multi-physics calculations of forming clusters, 
which include all four mechanisms for multiple formation. 


32. Disk Fragmentation 


3.2.1. Summary of the instability 


Independent of multiple formation on the scales of cores 
and filaments, protostellar disks, which form around the in- 
dividual stars, may also be prone to fragmentation into bi- 
nary or higher order multiples. The notion that gravitational 
instability in a shearing disk might produce binaries dates 
back over 30 years (Adams et al. 1989; Shu et al. 1990; 
Bonnell 1994). The physical mechanism at the heart of the 
instability is the same derived for the development of spi- 
ral arms in galaxies by Toomre (1964). The criterion for 
the growth of the axisymmetric, (m=0) mode of the insta- 
bility in razor-thin disks, where gravity overcomes support 
from thermal pressure on small scales and rotation on large 
scales is: 


(10) 


where c, is the disk sound speed, €? the disk orbital fre- 
quency, and X the surface density. The second form of the 
above equation relates the instability to the star-disk-mass 
ratio, M.. / M, and aspect ratio, H/r, where H is the disk 
scale height, and the order unity scale factor f depends on 
the assumed surface density power law. The Toomre Q cri- 
terion formally describes the onset of an instability rather 
than successful fragmentation of the disk into an indepen- 
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dent bound object. 

Numerical models of protostellar disk fragmentation of- 
ten reference a secondary criterion for instability to lead to 
fragmentation, the so-called “cooling criterion.” An early 
version of this mandate for disks to cool rapidly was first de- 
scribed by Thompson and Stevenson (1988) but was brought 
into much clearer focus by Gammie (2001). The constraint 
that disks radiate away heat on roughly the orbital timescale 
(teool «^3 — 71-1) is critical for understanding the transi- 
tion from the so-called gravito-turbulent state to fragmen- 
tation in viscously heated disks. For protostellar disks, the 
cooling constraint is often more easily satisfied than Q « 1 
(Kratter and Lodato 2016). Moreover for thick, massive 
disks, the critical Q value for fragmentation drops substan- 
tially below unity (Lau and Bertin 1978) due to the ex- 
tended vertical distribution of matter. Consequently, the 
canonical combination of Q~1 and teoo1 ~O(1) may be 
insufficient for predicting the onset of fragmentation in ir- 
radiated, massive disks (Tsukamoto et al. 2015; Takahashi 
et al. 2016). 


3.2.2. The onset of instability 


As the literature is replete with discussion of the nature 
of gravitational instability in general, we focus in this re- 
view on the likely pathway towards instability that is most 
relevant for multiple formation. Gravitational instability ( 
via modes m > 1) is most likely to arise in the outer disk, 
as Q(r)ocr?- (4*9)/?., where q is the power-law index of 
temperature T’xr~%, and p is the power-law index for sur- 
face density Har”. Except for pathological combinations 
of these indices, Q declines with radius. At disk radii of 
several 10s to 100s of au, the disk temperature is primar- 
ily set by stellar irradiation, and thus somewhat insensitive 
to changes in disk properties (Clarke 2009; Kratter et al. 
2008; Rice et al. 2011). Consequently, disks most likely 
become unstable through an increasing surface density, ©, 
that is not matched by an increased temperature. This sit- 
uation arises when the infall rate onto the disk exceeds the 
accretion rate through the disk and onto the star (Kratter 
et al. 2010; Harsono et al. 2011; Zhu et al. 2012). Such 
high infall rates are likely during the Class O phase. This 
scenario is particularly conducive to the formation of more 
equal mass binaries (rather than extreme mass-ratio BDs or 
super Jupiters) because the requirement for rapid infall also 
provides ample mass supply for the newborn fragment. In 
fact, for a wide range of disk conditions, the secondary ob- 
ject will out-compete the primary in accretion, driving mass 
ratios towards unity (Bate and Bonnell 1997; Ochi et al. 
2005; Young and Clarke 2015). This behavior is consis- 
tent with the preference for more equal mass ratios among 
close binaries (see $2.2). Note that even at arbitrarily high 
infall rates, the disk mass will not exceed the stellar mass, 
as global m —1 mode instabilities arise and lead to binary 
formation (Adams et al. 1989). Numerical simulations that 
invoke disk masses comparable to stellar masses as initial 
conditions are therefore somewhat unphysical. 
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An alternate pathway occurs when the disk temperature 
declines faster than accretion can alter the surface density. 
A likely candidate for this scenario is stellar luminosity 
variation due to accretion changes driven by gravitational 
or other instabilities. A sudden drop in the central luminos- 
ity could also precipitate fragmentation, so long as the drop 
in temperature does not correspond to a rise in the cool- 
ing timescale that prohibits fragmentation (Dunham et al. 
2014a; Kuffmeier et al. 2018). 


3.2.3. Fragment survival and subsequent evolution 


In order for gravitational instability and fragmenta- 
tion to produce binary formation, several benchmarks 
must be met. Fragments must cool efficiently so that 
they collapse to sizes well below their Hill radius, ry = 
a( Mirag/(3M.))'/3, before being sheared apart by inter- 
actions with the disk or even other fragments. Soon after 
formation, they are also subject to inward migration through 
the disk, which can lead to tidal disruption as the Hill radius 
shrinks faster than the object (Begelman and Cioffi 1989; 
Gammie 2001; Kratter and Murray-Clay 2011; Nayakshin 
2010; Boley et al. 2010). Merger with the host star is also 
possible (Vorobyov and Basu 2005a; Zhu et al. 2012). We 
defer a discussion of disk-driven migration in general to 
§3.3.2. Full modeling of these migration dependencies is 
crucial for predicting the expected population of binaries 
(separation, mass ratio) that derive from disk fragmentation 
plus migration. For example, their analytic model of disk 
fragmentation followed by inward migration and preferen- 
tial mass accretion onto the secondary, Tokovinin and Moe 
(2020) reproduced the observed excess twin fraction and 
BD desert among close solar-type binaries. 


3.3. Capture, Dynamical Evolution, and Migration 


The fragmentation modes discussed thus far provide ex- 
cellent explanations for a large swath of observed binary 
systems, even if they are not fully predictive analytic mod- 
els. One glaring omission in these models is the frequency 
and efficiency of orbital evolution post fragmentation. This 
omission is particularly acute for the closest binaries (P< 1 
au), where in-situ fragmentation is all but prohibited by the 
size of the first hydrostatic core (Larson 1969). To under- 
stand the origin of the closest binaries, as well as early evo- 
lution of orbital properties, we review the current models 
for capture and orbital evolution. 


3.3.1. N-body versus gas-mediated capture 


As described in §3.1, multiple systems can form at large 
distances, nearly independently as single stars, but never- 
theless be gravitationally bound to each other. Addition- 
ally, multiple stars in a star cluster may begin their lives un- 
bound, but due to energy and angular momentum exchange 
with the surrounding gas cloud or individual circumstellar 
disks, subsequently become bound and remain so. Other 
objects (often hierarchical multiples) will alternate between 
weakly bound and unbound over the first free fall time of 
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the cluster. We dub this process, seen in a range of scales 
in numerical simulations gas-mediated capture (Ostriker 
1994; Moeckel and Bally 2007; Bate 2012; Muñoz et al. 
2015; Geller et al. 2021b; Cournoyer-Cloutier et al. 2021). 
The net frequency of capture via these mechanisms is likely 
low; these systems would be indistinguishable observation- 
ally from those formed via fragmentation, except possibly 
cases where two disks appear to be on a collision course 
(see panel d, Fig. 5). 

On the other hand, in regions of high stellar density, in- 
teractions between truly unbound stars in the absence of 
gas can lead to the formation of new binaries and partner 
exchange. The outcome of such interactions is fairly well 
understood and can be modelled statistically. For exam- 
ple the adage that hard binaries are hardened by such in- 
teractions while soft binaries become softer (Heggie 1975) 
suggests that very frequent interactions should tend to drive 
binary separations towards a bimodal rather than unimodal 
distribution — notably inconsistent with most observations 
(see Fig. 2). More recent work has emphasized the im- 
portant contributions of higher order multiples to the over- 
all interaction rates in clusters of moderate density (Geller 
and Leigh 2015; Hamers and Samsing 2020). In partic- 
ular binary-single and binary-binary interactions can pro- 
duce a much more diverse set of outcomes, including the 
formation of very compact binaries (Dorval et al. 2017). 
While relevant for dense open clusters and globular clusters, 
these processes are sub-dominant in gas-rich environments. 
Moreover, even when rapid dynamical instabilities occur 
that change orbital parameters, system multiplicity, or hi- 
erarchy, they rarely lead to the order of magnitude changes 
in separation that are required to explain very close binaries 
(e.g., Wall et al. 2019). Orbital hardening is often driven by 
the ejection of a formerly bound companion; because the 
ejected object is typically the lowest mass, the total energy 
and angular momentum removed from the system in such 
an interaction causes more moderate changes in semi-major 
axis (Valtonen and Karttunen 2006; Kratter 2011). 

Partner exchange and shifting hierarchies also occur fre- 
quently in gas-rich environments (Bate 2012; Ryu et al. 
2017; Cournoyer-Cloutier et al. 2021). Such interactions 
are fundamental to shaping the final multiplicity distribu- 
tion observed for pre-MS and field stars and may even ex- 
plain some of the evolution between the Class 0 and I phase 
(see $4.2). Crucial as it may be, there is little fundamental 
theory to describe this process. Although some theories of 
multiple formation argue that most if not all objects are born 
in binaries (Kroupa 1995, 2008), modern numerical simu- 
lations produce some stars that are never bound to another 
cluster member (Bate 2012; Lee et al. 2019). 


3.3.2. Gas-driven migration 


In addition to allowing for capture, gas-rich star-forming 
environments also facilitate dramatic orbital evolution of bi- 
naries from their birth separation. The recent expansion of 
statistics for the youngest stars lends credence to the notion 
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that binaries form early and thus much of the orbital evolu- 
tion occurs before the natal cloud disperses (see $4.2). Mi- 
gration can occur either due to interactions with both bound 
or unbound low density molecular gas, or with a well de- 
fined circumbinary disk. 

Gas-dynamical friction: For widely separated bound bi- 
naries, the combined action of gas accretion and the gener- 
ation of a dense trailing wake provide a torque that reduces 
the semi-major axis. Numerous works have studied this so- 
called “dynamical-friction force" (Bate and Bonnell 1997; 
Ostriker et al. 1999; Stahler 2010; Lee and Stahler 2011). 
Lee et al. (2019) showed that a simple first order differential 
equation for à that includes both accretion and the torque 
due to dynamical friction with the gas =F x (—rhv,), can 
explain the observed rate of orbital decay in turbulent MHD 
simulations. The removal of angular momentum by mag- 
netic breaking also helps to drive inward migration (Zhao 
and Li 2013). Thus, so long as the gas cloud is present and 
stars continue to accrete, gas can efficiently reduce orbital 
separations from thousands of au to hundreds or even tens 
of au in well under a Myr. 

Though a simple model can explain the observed mi- 
gration in simulations, it does not easily translate into a 
generic prediction for the expected orbital evolution for a 
population of binaries, as it is a function of the local en- 
vironment. Nevertheless, phenomenological models have 
had some success in reproducing observed distributions 
(Tokovinin and Moe 2020). 

Disk-driven migration: Migration of binaries embedded 
in a shared circumbinary disk also crucially shapes the fi- 
nal orbital and mass ratio distributions. The classic theory 
for migration suggests that inward migration should dom- 
inate (Lubow and Artymowicz 1996), however this work 
only includes the torque from a circumbinary disk on point 
particles and neglects the often substantial torque generated 
by small individual disks surrounding each star, as well as 
the contribution from ongoing accretion. Accretion and cir- 
cumstellar torques in particular push objects in the oppo- 
site direction and can lead to substantial outward migration 
(Muñoz et al. 2019). The net torque including circumbinary 
disks, circumstellar disks, and ongoing accretion is a func- 
tion of binary mass ratio, eccentricity, effective viscosity, 
and aspect ratio (Satsuka et al. 2017; Dempsey et al. 2021). 
Ongoing work will map this space theoretically, which will 
inform future population synthesis models. However, it is 
clear that more extreme mass ratios experience inward mi- 
gration for typical disk properties. This is somewhat con- 
sistent with disk fragmentation simulations that undergo 
mergers, although numerical viscosity may also play a role 
(Vorobyov and Basu 2005b; Zhu et al. 2012). The produc- 
tion of near equal mass binaries at close separations may 
also derive from this process, as rapid inward migration at 
high mass ratios is followed by mass equalizing accretion 
from the circumbinary disk. This process can lead to stalled 
inward migration and even reverse the migration direction 
(Tokovinin and Moe 2020). 
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3.3.3. Secular evolution 


An alternative model for the origin of very close binaries 
(P <7 days) is Kozai-Lidov (KL) cycles coupled with tidal 
friction (Kozai 1962; Lidov 1962; Fabrycky and Tremaine 
2007). This model posits that binaries with P «10 days 
form from triple star systems with an inclined outer ter- 
tiary. Exchange of angular momentum between the or- 
bits can drive long timescale eccentricity-inclination oscil- 
lations. When the pericenter of the inner binary is pushed 
to very close separations, tides begin to circularize the orbit, 
both shutting off the oscillations and inducing drastic period 
decay (e.g., Dabringhausen et al. 2022). This model relies 
on the high fraction of triples amongst the closest binaries 
(Tokovinin et al. 2006, see $2.2). However, such oscilla- 
tions cannot account for the population of co-planar, close- 
in triples identified by Borkovits et al. (2016). A population 
synthesis study by Moe and Kratter (2018) found that less 
than half of binaries with P «10 days could be plausibly 
explained by KL cycles. One of the key reasons KL cycles 
do not reproduce the observed population is that the pre-MS 
close binary fraction is consistent with that of the field (see 
§2.4), suggesting that evolutionary mechanisms act on Myr 
not Gyr timescales: only a tiny fraction of KL susceptible 
systems have predicted orbital circularization timescales of 
less than a Gyr. The need for rapid migration likely requires 
gas-driven migration to efficiently drain energy and angular 
momentum from binary and triple star orbits. 


3.4. 


In past reviews, the limited number of star-cluster for- 
mation simulations and the small statistics of formed sys- 
tems therein have obscured trends in stellar multiplicity 
with birth environment. Over the past few years computa- 
tional advances have provided new insights into how molec- 
ular cloud properties affect stellar multiplicity by modelling 
more massive clouds, following larger dynamic ranges, and 
including more physical processes. While the statistics 
in such calculations remain modest, especially for very 
low (<0.1Mo) and high mass (> 10M.) stars, it is now 
clear that multiplicity exhibits more variation than the stel- 
lar IMF, which is surprisingly insensitive to environment 
(Offner et al. 2014). In this section, we discuss the relation- 
ship between environment, physical processes, and multi- 
plicity inferred from current state-of-the art star-cluster cal- 
culations. We present a summary of simulation properties 
and multiplicity statistics in Table 3. 


Impact of Environment on Multiplicity 


3.4.1. Radiation 


Radiative heating from protostars, which predominately 
influences gas within a few 10? au of the source in the 
case of low-mass stars, can significantly reduce disk frag- 
mentation ($3.2), thereby reducing the formation of brown 
dwarfs (Offner et al. 2009; Bate 2012). Bate (2019) mod- 
eled a 500M. unmagnetized, collapsing cloud with radia- 
tive transfer and a model for the diffuse ISM. The MF, sep- 
aration distributions, and mass ratios at the final time show 
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good agreement with the multiplicity of MS stars (e.g, see 
Fig. 6), despite the simulated stars being protostellar, i.e., 
significantly younger. Surprisingly, the same initial con- 
ditions evolved with a barotropic equation of state (Bate 
2012) produced statistically similar multiplicity with only 
a slightly higher MF for stars with M >0.2Mọ. The rela- 
tively small number statistics in the radiative transfer calcu- 
lation, combined with observational uncertainties, prohibit 
more rigorous conclusions about the impact of radiation on 
multiplicity. Note that the Bate (2019) calculations do not 
explicitly include heating due to stellar accretion and nu- 
clear processes, as the magnetized calculations from Cun- 
ningham et al. (2018) and Li et al. (2018) (Table 3), and 
therefore represent a lower limit on the impact of radiation 
on multiplicity. 

Models and simulations focusing on multiple formation 
via turbulent fragmentation, rather than disk fragmentation, 
find that radiative transfer has little impact on multiplicity. 
Core and filament fragmentation produce widely separated 
multiples, which form beyond the sphere of heating, such 
that calculations with and without protostellar heating give 
similar multiplicity statistics. However, even without disks, 
the statistics of very low-mass objects are over-predicted 
when heating is neglected, since turbulent fragmentation 
can proceed to smaller scales (Guszejnov et al. 2017). In 
simulations with disks, magnetic fields and thermal support 
due to protostellar heating both act to enhance disk stability 
and produce little small scale fragmentation (Offner et al. 
2016; Cunningham et al. 2018; Guszejnov et al. 2021). 


3.4.2. Protostellar outflows 


Protostellar outflows launched by accreting protostars 
entrain and expel dense core gas, both reducing the star 
formation efficiency of dense gas and injecting significant 
energy and momentum into the natal cloud environment 
(Bally 2016). Numerical simulations suggest that protostel- 
lar outflows can explain the offset between the stellar IMF 
and dense core mass function, helping to set the characteris- 
tic stellar mass (Offner and Chaban 2017; Guszejnov et al. 
2021; Mathew and Federrath 2021). Mathew and Federrath 
(2021) conducted a suite of MHD turbulent cloud simula- 
tions both with and without protostellar outflows. Runs in- 
cluding outflows reproduced the expected trend of increas- 
ing multiplicity with stellar mass but under-predicted the 
multiplicity of low-mass (M <0.03Mo) multiples. This 
deficit is likely due to the minimum calculation spatial reso- 
lution of ~ 10? au, which limits the formation of close sep- 
aration systems (although disk fragmentation does occur). 
Comparing the results of magnetized star cluster calcula- 
tions with and without outflows in Table 3 suggests outflows 
have little direct impact on CF and MF. 


3.4.3. Magnetic fields 


Magnetic fields provide additional pressure support that 
reduces gravitational collapse, thereby lowering the star for- 
mation rate, while increasing the characteristic stellar mass 
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TABLE 3 
STAR CLUSTER CALCULATION STATISTICS 


Properties? Statistics? 
Reference B R (0) C L (pc) M (Mo) Ax (au) MF CF à (au) dap (au) 

Bate (2019) "4 S 0.808 500 0.5 0.27-:0.03 — 0.503: 0.03 13 

Bate (2019) Z =0.01 Zo "4 S 0.808 500 0.5 0.29+0.05 0.52+0.04 9 

Bate (2019) Z=0.1Z6 v S 0.808 500 0.5 0.24 X 0.03 0.47 3: 0.04 8 

Bate (2019) Z =3 Zo ie v S 0.808 500 0.5 0.22+0.04 0.57 +0.03 25 tet 
Cunningham et al. (2018) p = 1.56 © v v v B 0.65 185 65 0.10 +0.03 0.18 +0.04 37 986 
Cunningham et al. (2018) p = 2.17 © v v v B 0.65 185 65 0.25+0.05 0.53 3: 0.09 126 564 
Cunningham et al. (2018) u = 23.1 * r4 v v B 0.65 185 32 0.20 + 0.06 0.38 +0.10 22 411 
Lee et al. (2019) p = 2 "4 B 2.0 601 50 0.14430.13 0.43+0.25 449 245 
Lee et al. (2019) yp — 8 v B 2.0 601 50 0.38+0.17 1.00+0.35 1146 1003 
Lee et al. (2019) w= 32 v das des B 2.0 601 50 0.08 +0.04 0.16 +0.06 1316 295 
Li et al. (2018) i, — 1.9 "4 "4 "4 B 4.2 3110 28 0.40+0.07 0.74+40.13 87 578 
Mathew and Federrath (2021) u — 6.6 v v B 2.0 TIS 100 0.342: 0.11 0.703: 0.25 291 560 


NOTE.—^ Physics included in the simulation: magnetic fields (B), radiative transfer (R), and outflows (O). Column C indicates a spherical (S) or box (B) cloud 
geometry; L, M, and Az indicate the cloud diameter, mass, and minimum gas spatial resolution, respectively. The metallicity, Z, is noted when it is not solar 
(Zo). Magnetized runs note the initial ratio of the gas mass to critical mass that can be supported by magnetic fields against collapse (mass-to-flux ratio"), ju. 
L= Meas/Merit = 2n V G Megas /Bz L?, where B, is the initial uniform magnetic field and L is domain length. ^ All entries except those of Mathew and Federrath 
(2021) assume binomial statistics and Poisson statistics to estimate the MF and CF uncertainties, respectively. Uncertainties for the Mathew and Federrath (2021) values 
adopt the standard deviation of the 10 runs. Bate (2019) uses closest separation to rank pairs and define stellar systems, while the other results use lowest binding energy 
to define bound systems. This choice has little impact on MF and CF, which differ by less than 1c. The last two columns list the median semi-major axis, à, and median 
2D pair separation computed by assigning system multiplicity following Tobin et al. (2022). © Statistics represent the final stellar distribution from a run with driven 
turbulence and one with identical initial conditions in which turbulence naturally decays. 


(Padoan et al. 2014; Guszejnov et al. 2021). While mag- 
netic fields might therefore naively be expected to decrease 
stellar multiplicity, they may instead enhance multiplicity. 
Offner et al. (2016) carried out a suite of magnetized tur- 
bulent collapse calculations and found that because mag- 
netic fields introduce asymmetry, which inhibits spherical 
collapse, they enhance turbulent fragmentation. Cunning- 
ham et al. (2018) and Lee et al. (2019) followed magne- 
tized turbulent clouds with different mass-to-flux ratios and 
showed that stronger field calculations produce comparable 
or higher MFs and CFs (see Table 3). Likewise, the MF and 
CF reported by Li et al. (2018) and Mathew and Federrath 
(2021) are higher than those of Bate (2019), which do not 
include magnetic fields. While strong magnetic fields re- 
duce protostar formation overall, they also produce fewer 
initial high-order systems. Consequently, fewer dynami- 
cal interaction occur, which would otherwise eject mem- 
bers and increase the fraction of single stars. This suggests 
a more subtle relationship between magnetic fields, stellar 
density and multiplicity. 

Magnetic fields treated in the ideal limit of perfect 
dust-gas coupling transport angular momentum efficiently, 
which can lead to compact (« 50 au) accretion disks or 
suppress disk formation altogether (Zhao et al. 2020). Spa- 
tial resolution of < a few au is required for disk formation 
in turbulent, ideal MHD calculations. Since no magne- 
tized cluster calculations listed in Table 3 adequately re- 
solve disks, multiples in these calculations formed primar- 
ily through turbulent fragmentation. 

Close binaries are under-represented in the resulting sep- 
aration distributions, which hints that disk fragmentation, 
in addition to turbulent fragmentation and dynamical in- 
teractions, is required to produce the observed population 
of close-separation field star binaries. However, Figure 6 
shows the resulting separation distributions of stars formed 
in some MHD calculations are consistent with the observed 
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distribution of wide-separation protostars from Tobin et al. 
(2022). This comparison is apt because the ages of the sim- 
ulated stars in these calculations are less than 1 Myr and, 
like protostars, they are less dynamically evolved than the 
population of field stars. Note that both simulations the pro- 
tostellar observations are incomplete for separations below 
7-25 au (see $4.2). Meanwhile the separation distribution 
from Bate (2019) is similar to the observed field population 
of solar-type stars and inconsistent with the protostellar dis- 
tribution despite simulated stellar ages 50.2 Myr. 

Non-ideal magnetic processes are important in regulat- 
ing angular momentum transport and setting disk proper- 
ties (Zhao et al. 2020, see also $3.2), however, most calcu- 
lations of star-cluster formation have assumed ideal MHD. 
In one exception, Wurster et al. (2019) modeled a suite of 
50M, magnetized turbulent clouds, including ambipolar 
diffusion, Ohmic dissipation, and the Hall effect. They find 
that non-ideal magnetic processes do not inhibit the forma- 
tion of multiple systems, which form independent of the 
initial field strength. Nevertheless, small number statistics 
(Nsys < 20) prohibit stronger conclusions about the impact 
of non-ideal effects on multiplicity. 


3.4.4. Turbulence 


Turbulent properties, such as the Mach number, virial 
parameter, and type of forcing (e.g., compressive versus 
solenoidal), influence the gas distribution and star forma- 
tion rate (Padoan et al. 2014; Federrath 2015). While 
the stellar IMF appears to be relatively insensitive to these 
details (Offner et al. 2014; Guszejnov et al. 2021), the 
impact of turbulence on stellar multiplicity is less clear, 
since few simulations have explored the impact of vary- 
ing turbulence properties on multiplicity. Cunningham 
et al. (2018) compared multiplicity in turbulent, magnetized 
clouds and found those with continuous turbulent driving 
and weaker fields (u= 1.56,2.17) over-predicted the mul- 
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tiplicity of solar-type stars, albeit the statistics were small. 
Mignon-Risse et al. (2021) found that the collapse of mas- 
sive cores with sub-Alfvénic turbulence (M Aequivv/v 4 < 
1), i.e., relatively stronger fields, produced single high-mass 
systems, while cores with super-Alfvénic turbulence pro- 
duced binary formation via disk fragmentation, which was 
driven by higher levels of angular momentum. 

A separate class of semi-analytical models focuses on 
the fundamental interaction between turbulence and grav- 
ity; this approach is simpler and computational cheaper 
than performing MHD simulations but cannot capture time 
dependence or the complex spatial relationship between 
turbulence and density. Despite a variety of approxima- 
tions, the Guszejnov et al. (2017) turbulence model (see 
$3.1.1) successfully predicts increasing multiplicity with 
stellar mass and the relative frequency of mass-ratios for 
both solar-mass and low-mass stars. Fig. 6 shows that the 
predicted MF agrees reasonably well with that of field stars, 
although the model over-predicts the MF of solar-type stars. 
This suggests that supersonic turbulent properties, such as 
the lognormal density distribution and power spectrum, to- 
gether with gravity can explain a variety of observed multi- 
plicity properties. However, Guszejnov et al. (2017) predict 
that the initial separation distribution peaks at ~ 600 au and 
has a deficit of close, d 10 au, binaries. This underscores 
that fully matching observations requires treatment of disk 
fragmentation and/or dynamical evolution (see $3.2). 


3.4.5. Metallicity 


It has long been understood that the evolution of metal- 
licity from the first Population III stars to today could sub- 
stantially impact the IMF (Abel et al. 2002). Gas metallicity 
influences the efficiency of cooling. Reducing the metal- 
licity of the low density, optically thin gas characteristic 
of star forming cores inhibits cooling via metal lines. On 
the contrary, for solar metallicity gas that is optically thick, 
as in protostellar disks, lowering metallicity can promote 
cooling by reducing the opacity. Bate (2019) carried out 
four calculations of forming star clusters with metallicity, 
Z, ranging from 0.01 to 3 times solar. The calculations in- 
cluded separate treatment for the gas and dust temperatures, 
but neglected radiative feedback and magnetic fields. While 
the shape of the stellar IMFs was largely independent of 
metallicity, the lower metallicity calculations experienced 
more fragmentation in cores, filaments, and disks due to the 
higher cooling rate of dense gas. This produced more close 
binaries (a S 10au), which varied from CF=0.1 +0.03 for 
3Zo to 0.53+0.12 for 0.01Zo. The close binaries in the 
lowest metallicity calculation also exhibited a small prefer- 
ence for lower mass ratios, which is likely caused by the en- 
hanced small scale fragmentation. However, Table 3 shows 
that the overall MF and CF values for the four metallicities 
fall within 1c of one another. 

As noted in $2.2.7, observations support a drastic de- 
crease in the close binary fraction as a function of metal- 
licity, with the binary fraction beyond ~ 200 au remaining 
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essentially unchanged (Moe et al. 2019; El-Badry and Rix 
2019). Analytic calculations suggest that the increase in 
only the close binary fraction is consistent with enhanced 
disk fragmentation, with wide binary formation mecha- 
nisms remaining essentially unchanged (Moe et al. 2019). 


4. Observations of Forming Multiple Systems 


Observations of multiplicity in the youngest stellar sys- 
tems are key to understanding their origins. Unlike the pre- 
MS and MS systems reviewed in $2.2, they have had little 
time to dynamically change, and therefore their properties 
more directly reflect the formation conditions. 

We define the youngest star systems, hereafter proto- 
stars, as stars that are still deeply embedded within a dense 
envelope of dust and gas that obscures them from direct ob- 
servation (i.e., Class O and I sources). Due to high opti- 
cal depths from the surrounding material, protostellar mul- 
tiplicity is best observed via high resolution interferometric 
observations at (sub)millimeter to radio wavelengths where 
the surrounding envelope is both optically thin and filtered 
out such that the circumstellar environment can be detected. 
At these wavelengths, young stars are primarily detected 
via thermal dust emission from circumstellar disks , ther- 
mal free-free emission from ionizing shocks at the base of 
jets. or from a combination of these processes (e.g., Segura- 
Cox et al. 2018; Tychoniec et al. 2018). 

In the Protostars & Planets VI review, discussion of pro- 
tostellar multiplicity centered on infrared observations of 
Class I systems, since there were few observations of mul- 
tiplicity in Class 0 systems (Reipurth et al. 2014). These 
studies tended to have inhomogenous separation distribu- 
tions, owing to different region distances, and many were 
limited to separations > 100 au, which does not probe the 
field star separation peak (see Fig. 6). Nevertheless, obser- 
vations provided tantalizing evidence that the multiplicity 
fraction of protostars is higher than that of pre-MS and MS 
stars (Looney et al. 2000; Chen et al. 2013) suggesting that 
most systems form as multiples. 

In this section, we review recent impactful case studies 
and large surveys using (sub)millimeter continuum obser- 
vations from interferometers that probe complete popula- 
tions of both Class 0 and Class I sources in nearby clouds. 
We describe new emerging trends and discuss their impli- 
cations for multiple star formation and evolution. 


4.1. 


Multiplicity can begin prior to the onset of star forma- 
tion when starless cores and filaments fragment to form 
multiple density peaks (see $3.1). We define starless cores 
as those with no detected protostar and independent of 
whether they are gravitationally bound. Observations of the 
internal structure of collapsing starless cores are required to 
understand how and when this mode of multiple formation 
occurs. Detecting multiple density peaks or extended inter- 
nal structure within a given starless core does not guaran- 


Starless Cores 
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axis distributions for simulated protostars.The number in each bin is normalized by the total number of singles and separations in each 
dataset. The “Li ea 2018 2D sep” curve shows the projected 2D separations of systems identified using the method of Tobin et al. (2022, 
see Table 3). Right: Multiplicity fraction versus primary system mass. Shading indicates bin ranges and binomial uncertainties. Black 


points show the bias-corrected MF for field stars (see Table 1). 


tee multiple formation, however, since the objects must also 
be self-gravitating and undergoing collapse. Current litera- 
ture tends to define any compact emission that is detected 
within a starless core using interferometry as substructure. 
The presence of two or more substructures suggests that the 
core is undergoing the first phase of collapse with possible 
incipient multiple formation, which is our main focus here. 
At the time of Protostars & Planets VI, surveys of star- 
less cores had found few cases of substructure (Schnee 
et al. 2010, 2012; Lee et al. 2013). The profusion of non- 
detections placed upper limits on substructure masses of 
~0.1 Mo. Deeper studies of a few objects successfully de- 
tected substructure: L 183 (Kirk et al. 2009), R CrA SMMI 
A (Chen and Arce 2010), Oph A SM1 (which has since 
been identified as protostellar, see §4.3.1), and Oph B2-N5 
(Nakamura et al. 2012). These deeper studies estimated 
substructure masses of ~0.001—0.01 Mo on sub-1000 au 
scales, in agreement with the earlier non-detections. 
ALMA has enabled large surveys of starless cores on 
sub-1000 au scales and down to the mass sensitivities of 
~0.001 Mo that are necessary to detect substructure (Dun- 
ham et al. 2016; Kirk et al. 2017). Dunham et al. (2016) 
found no substructure in 56 starless cores in Chameleon, 
whereas Kirk et al. (2017) detected substructure in only two 
out of 37 starless cores in Ophiuchus. The substructures de- 
tected by Kirk et al. (2017) include multiple density peaks 
in Oph A SMIN (see Fig. 5 and Friesen et al. 2018) and a 
single substructure in L1689N, a starless core to the east of 
IRAS 16293-2422 region (Lis et al. 2016). Oph B2-N5 was 
not detected by Kirk et al. (2017), however, likely due to a 
combination of lower sensitivities at long wavelengths and 
loss of emission from spatial filtering compared to Naka- 
mura et al. (2012). In a targeted study, Caselli et al. (2019) 
detected multiple density peaks in the starless core L1544. 
The relatively low substructure detection rate in starless 
cores implies that the cores have not reached sufficiently 
high central densities (possibly because collapse is not oc- 
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curring) or that the timescale over which the substructure 
is detectable before protostar formation is short. Kirk et al. 
(2017) reasoned that if the typical starless core lifetime is 
comparable to the free-fall time then two detections in 37 
starless cores implies a substructure timescale on order of 
5 kyr assuming a central density of 10° cm~?. Dunham 
et al. (2016) used synthetic observations of MHD simula- 
tions of turbulent cores and Bonnor-Ebert sphere models 
to constrain the likelihood of detection. They found that 
substructure is only detectable at their survey resolution to- 
ward the end of the starless stage (after > 0.1 Myr) when the 
central density exceeds ~ 107 cm-?, i.e., when the free-fall 
time itself is less than 10 kyr. 

Assuming typical core properties and that cores frag- 
ment on their free-fall timescale, Dunham et al. (2016) ex- 
pected roughly two cores in the Chameleon sample to show 
substructure. Most of these cores, however, are gravitation- 
ally unbound, suggesting they are less evolved and it is pos- 
sible some of them will not collapse to form stars at all. 
By contrast, the cores in Ophiuchus are denser, likely more 
evolved and therefore are more likely to contain substruc- 
ture (Kirk et al. 2017). Given the short timescales during 
which high densities are present prior to protostar formation 
and the large fraction of unbound cores in star-forming re- 
gions, detecting substructure in starless cores is challenging 
even with ALMA. Future studies will require large samples 
of cores or targeted observations with a balance of high res- 
olution and moderate scales. Offner et al. (2012b) found 
that synthetic observations of core fragmentation simula- 
tions that combined the main ALMA array with the ACA 
were key to detecting substructure in starless cores (see 
also, Dunham et al. 2016). 

It is worth noting that of the substructure detections in 
starless cores, most (4/6) show evidence of multiple density 
peaks in dust continuum, e.g., as shown in Fig. 5. This oc- 
currence rate is consistent with that of turbulent fragmen- 
tation models, which predict that a significant fraction of 
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cores and filaments that do collapse to form stars undergo 
fragmentation. This sample suggests that when substructure 
is detected there is a high probability that it shows evidence 
of multiple formation. 

The amount of core turbulence may be an important 
signpost for targeted studies searching for fragmentation. 
However, there has been no rigorous study comparing the 
turbulent properties of structured and unstructured cores. 
Figure 7 compares the non-thermal velocity dispersion, 
ont, to the virial parameter (a= Myi,/M) of starless 
objects that have detected substructure (black triangles) 
with the larger starless core populations in Ophiuchus and 
Chameleon (grey symbols). To first order, a core has super- 
sonic turbulence if ow /c, 7» 1 and subsonic turbulence if 
ONT/ 6s <1; it is bound if a= My, /M <2 and unbound if 
a= My /M »2 (McKee 1999). We calculate a virial mass 
assuming Myi, —5o, R/G, where oy is the total velocity 
dispersion, R is the size of the core, and we calculate c, 
assuming u= 2.37 and the gas is isothermal. We note that 
all values are taken from the literature where the turbulent 
properties, temperature, and mass estimates are all calcu- 
lated differently and thus a thorough consistent study of 
starless core properties is still needed. 
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Fig. 7.—: Starless core properties for L1688 in Ophiuchus 


(circles) and Chameleon (squares) with known structured starless 
cores as triangles. Data are primarily from single-dish telescopes 
and are calculated using different tracers, sensitivities, and tech- 
niques. We assume the core gas is isothermal to calculate the 
mass (measured from dust continuum data) and sound speed. The 
L1688 core properties are from Kerr et al. (2019) and use NH3 
observations. The Chameleon core properties use the N2H * ob- 
servations from Tsitali et al. (2015), with the temperature scaled 
to 12 K to match the dust masses from Belloche et al. (2011). The 
structured cores properties use data from the following references: 
L183 (Lattanzi et al. 2020; Karoly et al. 2020), R CrA SMM1 A 
(Harju et al. 1993; Bresnahan et al. 2018), Oph A SMIN (Pat- 
tle et al. 2015; Kerr et al. 2019), Oph B2-N5 (Motte et al. 1998; 
Friesen et al. 2009), L1689N (Pattle et al. 2015; Pan et al. 2017), 
and L1544 (Chacón-Tanarro et al. 2019; Koumpia et al. 2020). 


Figure 7 shows that the general starless core populations 
in Ophiuchus and Chameleon exhibit wide ranges of ve- 
locity dispersion and degree of gravitational boundedness, 
whereas the few cores that have detected internal structures 
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are all bound and most (5/6) have supersonic turbulence. 
The only subsonic core is L1544, which shows evidence 
of infall (Keto et al. 2015), indicating that it is collapsing. 
While cores without detections may still have structure be- 
low the observed sensitivity limits, the structured cores have 
density peaks that are more substantial and more likely to 
result in multiple star formation (e.g., Dunham et al. 2016). 
This distinction suggests that bound, turbulent, and collaps- 
ing cores are the best candidates for substructure searches. 


4.2. Protostellar Multiplicity Statistics 


Prior to Protostars & Planets VI the lack of large, un- 
biased protostellar samples led to overly broad conclusions 
based on small samples (Looney et al. 2000; Maury et al. 
2010). One of the first attempts to create a large sample was 
by Chen et al. (2013), with an archival study of submillime- 
ter observations toward protostars. This greatly increased 
the number of characterized systems, but the archival data 
had non-uniform sensitivity and resolution. 

The advent of ALMA and the upgraded VLA have now 
enabled truly uniform and unbiased surveys of multiplicity 
in nearby star-forming regions, where global statistics can 
be determined with high accuracy. Thus far, the protostel- 
lar populations of Perseus (Tobin et al. 2016b), Ophiuchus 
(Encalada et al. 2021), and Orion (Tobin et al. 2020, 2022) 
have been observed at uniform sensitivities and resolutions. 
This section reviews the statistical results of these censuses. 


4.2.1. Multiplicity and companion frequencies 


The MF and CF are key metrics for comparison to more 
evolved populations, but for a given sample they depend 
on separation range and protostellar Class, as well as any 
sample bias. For instance, the archival study of Chen et al. 
(2013) found MF=0.64+0.08 and CF=0.91-£0.05 for proto- 
stars with a separation range between 50 to 5000 au. How- 
ever, these values likely missed companions since the typi- 
cal spatial resolution was 100s of au rather than 50 au. 

Recent surveys of Perseus and Orion have well-defined 
angular resolutions and sensitivities (Tobin et al. 2022). 
From 20 to 10000 au, the MFs and CFs for the full pro- 
tostellar samples (Class 0, Class I and Flat Spectrum) 
are 0.30+0.03 and 0.44+0.03, respectively, in Orion and 
0.36+0.06 and 0.52+0.06, respectively, in Perseus. Both 
results are consistent within their uncertainties and both 
surveys have comparable minimum separations set by the 
angular resolution of the observations. For Class 0 proto- 
stars with separations between 20 to 10000 au the MFs and 
CFs are 0.38+0.05 and 0.62--0.05 for Orion and 0.47+0.09 
and 0.74+0.08 for Perseus, which are consistent despite 
environmental differences. These results are lower over- 
all with respect to Chen et al. (2013), despite the fact that 
Tobin et al. (2022) considers a wider range of separations. 
While the VANDAM surveys discovered a number of new 
companions with separations «200 au and wider compan- 
ions at 25,000 au, these were offset by the detection of 
many additional single systems, thereby reducing the total 
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MF and CF. Previous archival studies were limited to data 
biased toward the brightest millimeter sources, which also 
tend to have higher luminosities; Tobin et al. (2022) found 
that multiple protostar systems, regardless of separation, 
had systematically higher luminosities than single systems. 

Multiplicity statistics change between protostellar classes 
within the same 20 to 10^ au range of separations (Tobin 
et al. 2016b, 2022). In Orion, the MFs and CFs of Class 
I protostars are 0.23--0.04 and 0.32 + 0.05, respectively, 
while for Perseus they are 0.27-£0.09 0.35+0.09, respec- 
tively. Thus, Class I protostars have lower multiplicity than 
Class 0 protostars for the same separation range. The rel- 
ative MFs and CFs are ~40% and ~50% lower for both 
regions, respectively, having > 3oc differences. 

The MF and CF for Flat Spectrum protostars in Orion are 
consistent with the Class I statistics. Flat Spectrum proto- 
stars were not classified separately for Perseus but are com- 
bined with the Class I sample. These differences indicate 
multiplicity evolution between Classes (see $4.4.2). 

Other comprehensive studies of multiplicity have been 
carried out in Orion at infrared wavelengths, but these stud- 
ies were limited to «1000 au separations due to contamina- 
tion on larger scales. Kounkel et al. (2016) measured a CF 
of 0.14+9:011, which includes mostly Class I and Flat Spec- 
trum sources, demonstrating the large variation in MF and 
CF for different separation ranges and evolutionary stages 
(see §4.2.2). In this range of separations for a compara- 
ble sample of Class I and Flat Spectrum protostars, Orion 
has a MF and CF of 0. 100.02 and 0.1 10.02 (Tobin et al. 
2022), respectively. These are slightly lower but still con- 
sistent with the statistics from Kounkel et al. (2016), which 
is likely because the centimeter/millimeter observations are 
more incomplete toward Class I/Flat Spectrum protostars 
than near-infrared observations (Tobin et al. 2022). 

Finally, Encalada et al. (2021) characterized the Class 
I and Flat Spectrum protostellar multiplicity in Ophiuchus 
for the separation range between 15 au and 1600 au. They 
found a combined Class I and Flat Spectrum MF and CF of 
0.25+0.09 and CF=0.33+0.10, respectively, which is con- 
sistent with those of Orion and Perseus. 


4.2.2. Separation distributions 


High spatial resolution, uniform surveys now measure 
companion separation distributions down to <50 au, simi- 
lar to studies done for MS stars (see $2). The large sam- 
ple sizes enable studies as a function of protostellar class, 
which have distributions that appear distinct from those of 
MS stars (Figures 2 & 6). 

Tobin et al. (2016b) detected an apparent bimodal sepa- 
ration distribution for the full sample of Perseus protostars 
with peaks at ~75 au and ~3000 au. The positions of the 
two peaks indicate both disk and turbulent fragmentation 
operate. However, the wide-separation feature is mainly 
produced by the Class O sources, and a reanalysis of the 
Perseus data does not rule-out a log-flat distribution. Orion 
also exhibits a bimodal separation distribution with similar 
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peak locations (Tobin et al. 2022). The larger samples for 
Orion alone and the combined Orion and Perseus samples 
strongly disfavor a log-flat distribution for both Class 0 and 
Class I protostars. 

(Encalada et al. 2021) measured the separation distribu- 
tion from 15 to 1600 au for Class I and Flat spectrum proto- 
stars in Ophiuchus. Despite small numbers, this distribution 
is quite similar to that of Perseus and Orion and has an aver- 
age separation of 183 au. The similarity of separation dis- 
tributions in different star-forming environments suggests 
multiple systems form via similar physical processes. 


4.3. Markers of Binary Formation 


4.3.1. Multi-generational systems 


MS stars in multiple systems are generally expected to 
have formed together and have similar ages. Nevertheless, 
theoretical models predict age differences up to a few Myr 
between stars in multiple systems. While a few Myr age 
difference is difficult to identify in MS multiples, it can po- 
tentially be observed in younger systems. In this section, 
we review forming multiple systems that appear to have 
distinct star formation epochs and discuss the degree of age 
difference. We use the term non-coeval to describe apparent 
protostellar age differences. 

Murillo et al. (2016) examined the spectral energy distri- 
butions (SEDs) of 16 protostellar multiples in Perseus with 
data from micron to millimeter wavelengths and found that 
roughly a third of the wider (= 2000 au) multiple systems 
contained members with different Classes. Assuming Class 
is an accurate proxy for age, the implied age differences 
are consistent with those predicted (At <3 x 10° Myr) by 
numerical simulations (e.g., Lee et al. 2019). 

Recent ALMA observations have identified starless sub- 
structures near protostars that may point to core fragmen- 
tation after star formation is well underway. Two instances 
are Oph A SM1 and Oph A N6 (Friesen et al. 2018). Oph 
A N6 exhibits two peaks, one of which appears protostel- 
lar (Friesen et al. 2018). Recent ALMA observations show 
that one of the peaks in Oph A SM1 is also likely protostel- 
lar (Friesen et al. 2014). Oph SM1 appears connected to the 
structured starless core SMIN to the north (see Fig. 5) and 
a chain of other protostars to the south (Friesen et al. 2014; 
Kirk et al. 2017), suggesting that the entire region may be 
undergoing filamentary fragmentation with an age gradient 
from north to south. 

Another example is B5, a bound elongated core that con- 
tains two starless filaments that extend away from a central 
Class I protostar (Pineda et al. 2011). Pineda et al. (2015) 
found three dense condensations within these filaments as 
shown in Fig. 5, which are incipient sites of star formation 
that could ultimately produce a non-coeval quadruple sys- 
tem. Schmiedeke et al. (2021) showed that the mass per unit 
length (M/ L) exceeds the expected critical value required 
for stability (see 3.1), in agreement with the presence of 
the embedded condensations. Thus, core and filament gas 
surrounding young protostars may still fragment, producing 
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age differences comparable to the duration of the protostel- 
lar lifetime, i.e., ~0.5 Myr. 


4.3.2. Protostellar disk stability 


Observations of massive circumbinary and circum- 
multiple disks provide evidence for the occurrence of disk 
fragmentation. A number of studies have found extended 
structures around already formed, tight multiple systems 
(e.g., Tobin et al. 2016b,a; Lim et al. 2016; Brinch et al. 
2016; Takakuwa et al. 2017; Kraus et al. 2017; Harris et al. 
2018; Artur de la Villarmois et al. 2018; Alves et al. 2019; 
Lee et al. 2020a; Tobin et al. 2020), some of which exhibit 
Keplerian rotation signifying that they are circumbinary 
disks (e.g., Murillo et al. 2016; Artur de la Villarmois et al. 
2019; Reynolds et al. 2021). Direct confirmation of disk 
fragmentation, however, is difficult since protostars formed 
by turbulent fragmentation combined with orbital evolution 
can drive wide companions to closer separations (see $3.3). 
The Toomre Q parameter provides insight into the actual 
degree of disk stability ($3.2). In this section, we adopt a 
critical Toomre Q of Qerit = 1; while this threshold for frag- 
mentation is not universal, the variation is far smaller than 
the observational error bars (Kratter and Lodato 2016). 

The number of disk stability measurements using the 
Toomre Q parameter are increasing, although all calcula- 
tions of Toomre Q should be considered first-order esti- 
mates due to uncertain temperatures, dust opacities, and 
dust-to-gas conversion factors that are necessary to obtain 
surface densities. Additional diagnostics, such as evidence 
of disk structure, are therefore necessary to verify instabil- 
ity. The circum-multiple disk of L1448 IRS3B in Perseus 
(Fig. 5) is a prime candidate for disk fragmentation. The 
disk geometry and rotation is centered on a close binary 
system, but a third object is detected in the outer spiral arm 
(« 300 au from the binary). This tertiary feature also ap- 
pears protostellar, with a compact source that is driving an 
outflow (Reynolds et al. 2021). At the radius of the tertiary 
source the disk has Q « 1 (Tobin et al. 2016a; Reynolds et al. 
2021), supporting the theory that this protostar companion 
was produced by disk fragmentation. Similarly, the Class 
I protostar HH111 MMS has spiral structure, Q «1 (Lee 
et al. 2020a) and a wide companion. Although no fragments 
are detected toward this disk yet, the similarities to L1448 
IRS3B suggest it is a prime candidate for fragmentation. 

Nevertheless, most protostellar disks with measured 
Toomre Q values appear to be stable. Tobin et al. (2020) 
measured the Toomre Q parameter for 259 disks in Orion 
that are in either single star systems or in wide binaries, 
finding that only six were consistent with being gravita- 
tionally unstable (see also, Sharma et al. 2020). Moreover, 
some structured circumbinary disks have large Toomre Q 
values indicating that they are stable. Reynolds et al. (2021) 
found hints of spiral structure toward the disk of L1448 
IRS3A, a more distant companion to L1448 IRS3B (Tobin 
et al. 20162), with Qz5. Similarly, Alves et al. (2019) 
and Diaz-Rodriguez et al. (2022) each found prominent 
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spiral structure in a Class I circumbinary disk and Q = 4. 
These cases show that spiral structure alone does not indi- 
cate gravitational instability. Alternatively, spiral structure 
may be caused by dynamical processes such as interactions 
between the embedded protostars or it may be a remnant 
of a previously more massive disk from which the stars 
formed. We note that protostellar disks can accrete material 
from the surrounding envelope. Thus, this simple picture of 
gravitational instability is more complicated, since the disk 
mass and surface density change over time. 

The lack of marginally-unstable disks suggests disk in- 
stability is either short-lived or rare. In the former case, 
disks fragment quickly such leaving behind a compact mul- 
tiple system. In the latter case, if few disks ever become 
gravitationally unstable then most protostellar companions 
at separations <300 au must be produced by migration. 
Further observations tracking fragmentation and disk sta- 
bility in different cloud environments across multiple evo- 
lutionary stages are necessary to distinguish the fractions of 
close multiples produced in situ within disks versus those 
produced by migration. 


4.3.3. Protostellar disk properties 


In this section, we compare protostellar disk properties 
in systems with multiple and single stars. For simplicity, we 
exclude systems that have confused disks that are blended 
with circumbinary structures (see §5.2) or very close com- 
panions. Fig. 8 shows cumulative mass and size distribu- 
tions for protostellar disks in Orion, which is the largest 
available sample (Jobin et al. 2020). The disks are sepa- 
rated into three populations: single-star systems, close and 
intermediate multiple systems, and wide multiple systems. 

Disks in close and intermediate multiple systems are 
statistically smaller in mass and size than those in wide 
multiples with p<0.01 from a Kolmogorov-Smirnoff test, 
whereas disks in single and wide-multiple systems tend to 
be consistent in both mass and radius. These trends indicate 
that the proximity of companions in close and intermediate 
multiple systems affects young disk properties. Wider com- 
panions, however, have little to no effect, and disks in those 
systems appear similar to those of single stars (see §5.1). 


4.3.4. Protostellar system orientations 


Stars that form via turbulent fragmentation should have 
uncorrelated angular momenta, since these protostars form 
in distinct, well-separated core substructures, whereas com- 
panions produced by disk fragmentation, form in a single 
plane and should have preferentially aligned angular mo- 
mentum vectors (Offner et al. 2016; Bate 2018). Early ob- 
servations of field star rotation appeared to support this pic- 
ture (Hale 1994). However, a rigorous statistical reanalysis 
of the Hale (1994) data concludes that these data are insuf- 
ficient to make any claims about spin-alignment (Justesen 
and Albrecht 2020). A handful of individual measurements 
of eclipsing binaries reveal heterogeneous orientations, but 
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Fig. 8.—: Comparison of circumstellar disk mass (left) and radius (right) for single and multiple Class 0 and I protostellar systems in 


Orion Tobin et al. (2020). Masses are calculated assuming optically thin dust emission at 870 um and an assumed dust temperature of 
20 K and disk radii are derived from the 2c extent of Gaussian fits to dust continuum images. The division at 300 au contrasts the effect 
of intermediate and wide separation companions on disk masses and extents. Surprisingly, systems with wide companions tend to have 
larger and more massive disks than those without companions. The maximum fraction is «1.0 primarily due to unresolved disks. 


with a preference towards alignment (Sybilski et al. 2018). 
To complicate matters, a myriad of mechanisms can torque 
multiple star systems into or out of alignment, such as tur- 
bulent accretion, magnetized accretion, ejections, inclina- 
tion oscillations, and damping from the disk (Lubow and 
Ogilvie 2000; Valtonen and Karttunen 2006; Batygin and 
Adams 2013; Fielding et al. 2015; Lee et al. 2017b; Jen- 
nings and Chiang 2021). 

Observations of forming systems provide a stronger test 
of the natal spin orientations of multiple systems before 
dynamical processes affect them. A variety of recent sur- 
veys have explored the angular momentum of young mul- 
tiple systems by using the outflow or disk orientations as a 
proxy for angular momentum direction. We note that even 
in these studies, where an angular momentum direction can 
be robustly measured and assigned to each protostar, such 
proxies only reflect the current angular moment direction of 
gas on ~ 1— 100 au radii from the stars and do not directly 
probe the (proto)stellar spin or system orbit. 

Using observations of outflows, Lee et al. (2016) studied 
the orientations of nine multiple systems in the Perseus 
cloud with separations 2000 au using data from the 
MASSES survey (Stephens et al. 2019). They found that 
the relative outflow orientations in these systems were con- 
sistent with either random or preferentially perpendicular 
alignments. For companions with separations less than a 
few hundred au, however, individual outflows can overlap 
in projection (e.g., Stephens et al. 2017). Blended outflows 
are difficult to disentangle and can bias the angular differ- 
ences toward misalignment (Offner et al. 2016). 

To date, there are no systematic studies of protostellar 
disk alignments. Many protostellar binaries with separa- 
tions less than a few hundred au have disks with similar 
inclinations and position angles (see Tobin et al. 2016b, 
2020). This similarity in disk geometries implies aligned 
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angular momentum vectors, in agreement with the simple 
picture of disk fragmentation. Similarly, several wider mul- 
tiples have misaligned disks in agreement with turbulent 
core fragmentation producing random angular momentum 
vectors, i.e., NGC 1333 IRAS2 (Tobin et al. 2015), IRAS 
04191+1523 (Lee et al. 2017a), BHR 71 (Tobin et al. 2019), 
GSS 30 (Sadavoy et al. 2019), HOPS 182 (Tobin et al. 
2020), and HH111 MMS (Lee et al. 20202). 

Nevertheless, some protobinary systems with separa- 
tions < 100 au seem misaligned. Oph-IRS 43 and Oph-IRS 
67 are close binaries that appear misaligned with respect to 
circumbinary disks, whereas the circumstellar disks for the 
tight binary in L1551 IRS 5 are inclined moderately (~ 30?) 
relative to the circumbinary disk and the circumstellar disks 
in a wider binary, L1551NE, are inclined relative to each 
other (Brinch et al. 2016; Takakuwa et al. 2017; Artur de 
la Villarmois et al. 2018; Cruz-Sáenz de Miera et al. 2019). 
Recent high resolution ALMA observations of the outflow 
of VLA 1623A show that its tight binary (< 30 au) has mis- 
aligned outflows and disks (Hara et al. 2021). Similarly, 
a disk resolved in the triple system IRAS 16293-2422 ex- 
hibits a ~ 90° misalignment with the surrounding circumbi- 
nary disk (Maureira et al. 2020). Such misaligned systems 
may have formed at larger separations via turbulent frag- 
mentation and then migrated inward (see $3) or the system 
may have been perturbed from an originally aligned orien- 
tation by dynamical interactions. 


4.4. Environment Around Forming Binaries 


4.4.1. Core and envelope properties 


Proto-multiple systems are embedded in dense cores and 
envelopes, which likely provide clues about multiple forma- 
tion. However, to date there has been few systematic studies 
of the surrounding environments of such systems. Sadavoy 
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and Stahler (2017) compared the separations of multiple 
protostellar systems in the Perseus cloud with the proper- 
ties of their host cores. They found that the projected dis- 
tance vectors of wide binaries with separations = 500 au are 
preferentially aligned with the long axes of their host cores, 
whereas tighter systems are more randomly orientated. This 
correlation suggests core fragmentation proceeds within a 
preferred plane, possibly due to rotation, magnetic fields, or 
both (see $3), whereas disk fragmentation or dynamics ran- 
domizes the orbit orientation relative to the host core. Luo 
et al. (2022) examined the binarity and core properties of 43 
systems in the Orion cloud. They found that the cores host- 
ing multiple systems tended to have higher column and vol- 
ume densities and higher Mach numbers than those hosting 
single stars. Higher densities and Mach numbers promote 
thermal Jeans fragmentation and turbulent fragmentation, 
respectively. This result suggests that multiple stars form 
from cores that are both turbulent and dense, in agreement 
with the structured detections for starless cores (see Fig. 7). 
Magnetic fields can greatly influence the formation of 
multiple systems, particularly if the field is aligned with the 
system angular momentum vector (see $3). However, mag- 
netic field observations exist for only a few protobinary sys- 
tem envelopes (< 1000 au). L1448 IRS 2 is a Class 0 proto- 
binary system with an inferred envelope magnetic field re- 
sembles an hourglass and is aligned with the system rotation 
axis (Kwon et al. 2019). Models indicate this field geometry 
enhances magnetic breaking and might produce more com- 
pact disks. Nevertheless, L1448 IRS 2 has a circumbinary 
disk (Tobin et al. 2018). Similar polarization morpholo- 
gies are seen toward the envelopes or circumbinary material 
of the protobinaries Per-emb-2 (Cox et al. 2018), BHBOT- 
11 (Alves et al. 2018) VLA 1623A (Sadavoy et al. 20182), 
IRAS 16293A (Sadavoy et al. 2018b), and BHR 71 (Myers 
et al. 2020), although it is unclear whether the polarization 
uniquely traces pinched fields versus outflow cavities. 
Some protobinaries also have complex circum-multiple 
material. Theoretical studies predict that bridges of dust and 
gas can connect binary stars, thereby affecting their forma- 
tion and accretion properties (83.1.2). At present, only the 
IRAS 16293-2422 protobinary system shows such a fila- 
ment of dust and gas connecting young low-mass stars (sep- 
arated by ~ 700 au Jørgensen et al. 2016). This dust bridge 
has relatively uniform polarization indicating a structured 
magnetic field (Sadavoy et al. 2018b), allowing gas to flow 
freely along the bridge. MHD simulations by Kuffmeier 
et al. (2019) show that such features occur when binaries 
form at wider separations (> 1000au) and migrate inward. 
If bridge features are transient, large high-resolution sur- 
veys are needed to constrain their properties and lifetimes. 
Binary interactions may also enhance the degree of bursty 
accretion and affect gas chemistry (Jørgensen et al. 2022). 


4.4.2. Stellar densities 


The local density of stars influences the evolution of 
multiple systems, such that wider members may be dynam- 
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ically stripped from systems located in dense cluster envi- 
ronments. This process is illustrated by the deficit of Class 
II and III companions with separations greater than 250 au 
in the ONC, where the stellar surface density is 770 pc? 
(Reipurth et al. 2007). The companion frequency is higher 
for field stars and pre-MS stars in Taurus, which have stellar 
densities down to 35 pc? (see right panel of Fig.2). 

Kounkel et al. (2016) completed a large Hubble Space 
Telescope and NASA Infrared Telescope Facility Survey 
of Class I and Flat Spectrum multiplicity toward proto- 
stars in the Herschel Orion Protostar Survey (HOPS, Furlan 
et al. 2016). They estimated the YSO density around each 
protostar by combining the HST data with complementary 
Spitzer and X-ray data. Adopting 45 pc? as the boundary 
between high and low YSO surface density, they find the 
CF for the 100-1000 au separation range is larger by ~50% 
(0.05 to 0.11) in higher surface density regions. 

Building on these results, Tobin et al. (2022) used the 
ALMA and VLA surveys of Orion to examine the influence 
of YSO surface density on protostellar multiplicity. They 
find higher MF and CFs in regions with YSO surface densi- 
ties above 30 pc ?, confirming the results of Kounkel et al. 
(2016). They also demonstrated a relationship between stel- 
lar density, separation range, and protostellar Class. For a 
separation range of 20 to 10^ au, all protostellar Classes 
have higher MFs and CFs in higher YSO density regions. 

This trend holds for Class I and Flat Spectrum sources 
with smaller separation ranges of 20 to 10? au and 100 to 
10? au. Class 0 protostars, however, have MFs and CFs 
that are statistically consistent in low and high YSO surface 
densities over all separation ranges. 

These results indicate that YSO surface density impacts 
multiplicity on scales « 1000 au, but that impact is only ap- 
parent after the Class 0 phase. Tobin et al. (2022) suggested 
that the Class O results reflect the primordial multiplicity, 
which does not have a strong dependence on YSO density. 
In contrast, the elevated MFs for more-evolved protostars 
are likely produced by the migration of companions from 
71000 au to «1000 au (see $3.3). The higher CFs from 
20 to 10^ au for all Classes at higher YSO density natu- 
rally provide a reservoir of companions that can migrate to 
smaller separations. This is consistent with the turbulent 
fragmentation model, which predicts companions form at 
wide separations and then migrate inward ($3.1). 

Thus, higher YSO densities may enhance rather than re- 
duce multiplicity. These high YSO densities are likely pro- 
duced by higher initial gas densities, which may favor frag- 
mentation (Gutermuth et al. 2011). We note, however, that 
the YSO surface densities in the Orion molecular clouds 
are lower than the stellar densities in the ONC, which ex- 
ceed 100 pe~?. Above some threshold, high YSO density 
likely hinders multiple star formation and the ability of such 
systems to retain wide companions. 
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5. IMPACT OF MULTIPLICITY ON PLANET FOR- 


MATION AND SYSTEM ARCHITECTURES 


The high frequency of multiple systems has the potential 
to have a strong impact on planet formation, both in terms 
of the overall frequency of planetary systems and on the ar- 
chitecture of those systems. In this section, we review the 
connections between multiplicity and planet formation, but 
we do not attempt a comprehensive review of all observa- 
tions of protoplanetary disks in multiple systems. Instead, 
we focus on some specific areas where the observations 
place particular constraints on planetary system properties, 
either current or future. After discussing disks in binaries, 
we summarize the final population statistics for planets in 
binary systems among field stars, again focusing on aspects 
that shed light on formation models. 


5.1. Raw Material for Planet Formation: Disk Sizes 
and Masses in Multiple Systems 
5.1.1. Disk masses and millimeter fluxes 


It has long been known that stellar companions with 
separations on the order of disk radii (~ 10 to 100 au) af- 
fect disk structure, reducing disk masses compared to disks 
around single stars (Beckwith et al. 1990; Osterloh and 
Beckwith 1995; Jensen et al. 1996; Harris et al. 2012). In- 
deed, surveys that start from a sample of young stars with- 
out pre-selecting them based on whether or not they host 
disks find that the overall disk frequency is lower in close 
binaries than in wider binaries or single stars (e.g., Cieza 
et al. 2009; Cheetham et al. 2015). 

Since Protostars & Planets VI, ALMA observations 
have added nuance to such studies due to the combination 
of ALMA’s high spatial resolution and sensitivity, particu- 
larly focusing on stars that still retain disks, i.e., the Class II 
phase, and measuring the properties of those disks that still 
survive at this age. ALMA has resolved many more binaries 
at millimeter wavelengths, allowing comparative studies of 
disk radii, and determining the flux distribution in circum- 
primary and circumsecondary disks. Increased sensitivity 
has allowed millimeter-wavelength detection of many more 
disks around secondary stars in particular, enabling compar- 
isons of disk properties around primary and secondary stars 
within binary systems. Such studies also extend to much 
lower stellar masses for both binary and single stars, thus 
including disks that are arguably more typical of the overall 
population of young stars. 

Akeson et al. (2019) combined new ALMA observations 
with literature data to assemble a sample of Class II bina- 
ries in Taurus wider than 30 au and with spectral types of 
M6 and earlier and compared these to a sample of single 
stars. To isolate the effects of binarity as cleanly as possi- 
ble, the single-star sample was required to have high spa- 
tial resolution in the optical or IR, and the binary sample 
was restricted to double systems only, removing the com- 
plicating effects of hierarchical systems. As has been seen 
previously, binaries with a companion closer than 140 au 
show reduced mm flux, while the mm flux distribution in 


25 


The Origin and Evolution of Multiple Star Systems 


wider binaries is indistinguishable from that of single stars. 
However, some interesting new results emerge when the in- 
dividual stellar masses are taken into account. It is well 
established that disk mass is correlated with stellar mass 
(e.g., Andrews et al. 2013; Pascucci et al. 2016). To account 
for this, Akeson et al. (2019) compared the distributions of 
Mais / M. in binaries and single stars and found that both 
primaries and secondaries, in both close and wide binaries, 
have reduced disk masses (on average) compared to sin- 
gle stars when corrected for stellar mass (Fig. 9a). While 
the difference between wide binaries and single stars is less 
than that for close binaries vs. singles (Fig. 9b), this result 
suggests that even the disks in wider binaries are influenced 
by membership in a binary system. In addition, after cor- 
recting for stellar mass, primaries do not dominate the disk 
mass in a given binary system; the ratio Ma; / M. has 
about the same distribution for primaries and secondaries. 
We note here that this and other similar studies assume that 
the mm continuum emission is largely optically thin, so that 
the emission traces the total dust mass in the disks. 

Zurlo et al. (2020, 2021) studied disks in binaries in o 
Oph and Lupus, combining new adaptive optics observa- 
tions to search for companions in their sample with previ- 
ous ALMA observations of disks from the ODISEA sur- 
vey (Cieza et al. 2019; Williams et al. 2019). They found 
that the most massive disks are present only around single 
stars, and they argue that for more typical disks (those with 
Maust < 90M3) the dust mass distributions are similar for 
single stars and multiple systems. 

While the studies discussed above focus on Class II sys- 
tems in 1-3 Myr old star-forming regions, Barenfeld et al. 
(2019) find somewhat different results for Class II systems 
with ages of ~5—11 Myr in Upper Sco. Unsurprisingly the 
0.88 mm flux density distributions of all groups (close bi- 
naries, wide binaries, and single stars) in Upper Sco are 
shifted to lower values compared to the younger Taurus 
association. In contrast to Taurus, Upper Sco exhibits no 
significant difference in the mm flux distributions of these 
three groups; in particular, close binaries have the same flux 
distribution as single stars. The stellar mass distributions of 
the three groups are similar as well, so this appears to be a 
real difference in the disk properties and not an artifact of 
different stellar masses. 


5.1.2. Disk radii 


One notable change in the ALMA era, compared to work 
with previous millimeter facilities, is the ability to measure 
(or at least constrain much better) disk radii, not only for the 
largest and brightest disks, but for more typical disks across 
arange of stellar masses. Cox et al. (2017), studying a sam- 
ple of mostly Class II stars in p Oph, found that dust disk 
radii among the binaries in their sample are smaller on av- 
erage than those of single stars. Manara et al. (2019) found 
a similar result in Taurus, with the added information that 
their single-star and binary samples have similar distribu- 
tions of stellar masses, showing that this result is not a side 
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Fig. 9.—: Disk mass (assuming optically thin emission) for Taurus Class II binaries and single stars normalized by stellar mass. Left: 


After normalizing by stellar mass, secondary and primary stars in binaries have similar disk mass distributions, but both are lower than 


the distribution of single star disk masses. Right: Smaller-separation binaries have lower-mass disks on average (relative to their stellar 


masses) than wider binaries, which in turn are lower than that of single stars. (Figure adapted from Akeson et al. 2019.) 


effect of the single stars in the sample being preferentially 
more massive. In contrast, Barenfeld et al. (2017) found 
that dust disks in Upper Sco in general, i.e., in both singles 
and binaries, are smaller by a factor of ~3 than those in 
similar-stellar-mass samples in Lupus or Taurus and Ophi- 
uchus (Tazzari et al. 2017; Tripathi et al. 2017, see also 
Hendler et al. 2020). Barenfeld et al. (2019) suggested that 
much of the difference between the single and multiple-star 
samples may be due to differing evolution of disk radii; tidal 
truncation may limit young disk sizes in binaries, but single 
star disks later “catch up” in reducing their sizes, perhaps 
by radial drift of dust, photoevaporation, dust fragmenta- 
tion, and/or viscous accretion (e.g., Gorti et al. 2015). 
Both Cox et al. (2017) and Manara et al. (2019) found 
that the measured dust disk sizes in binary systems are 
smaller than would be expected from tidal truncation alone, 
unless all systems had very large eccentricities. However, 
dust and gas disk sizes need not be the same; Zagaria et al. 
(2021b) showed that once the faster radial drift of dust ex- 
pected in binary disks is taken into account, the observed 
dust disk sizes in binaries are consistent with the expected 
tidal truncation radii for modest eccentricities. Rota et al. 
(2022) confirmed this observationally, showing that most of 
these systems have gas disks that are larger than their dust 
disks, consistent with radial migration of the dust. There 
is no discrepancy between the observed gas disk radii and 
expectations from tidal truncation models, assuming eccen- 
tricities consistent with the field eccentricity distribution. 
For the older disks in Upper Sco, it is not only curious 
why the disks in binaries have fluxes similar to those for 
single stars, but why (at least for the binaries) the disks are 
there at all. As noted above, at 1-2 Myr ages the binary 
disks were already small, and both gas dispersal and radial 
drift of the dust should be faster in binaries, with models 
predicting that disks in binaries with separations less than 
100 au should be depleted completely in just a few Myr 
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(Rosotti and Clarke 2018; Zagaria et al. 2021a). There 
is also no evidence of a circumbinary reservoir that is re- 
plenishing these disks (see §5.2). Here it is important to 
note that the studies of Barenfeld et al. (2017) and Baren- 
feld et al. (2019) (as well as most of the other studies cited 
above) adopt a target sample containing only Class II sys- 
tems in order to study why systems with ages similar to 
Upper Sco still retain some disk material. Thus, the simi- 
larity of the disk properties in binaries and singles may re- 
flect similar mechanisms for retaining disks on timescales 
of 10 Myr in any type of system. If the remaining disks 
are those that have developed large pressure bumps and/or 
dust traps, the similar disk properties might suggest that the 
dust-retention mechanism, rather than multiplicity, is now 
the dominant force driving disk evolution. 


5.2. Circumbinary Disks 


Most of the above discussion concerns disks that are cir- 
cumstellar with a given disk centered around one specific 
star in the system, but orbital dynamics also allow for cir- 
cumbinary disks, the presumed birthplace of circumbinary 
planets, roughly a dozen of which have now been discov- 
ered by Kepler (summarized in Martin 2019) and TESS 
(Kostov et al. 2020, 2021). 

While there are some notable examples of Class II cir- 
cumbinary disks, they are surprisingly rare. Stable cir- 
cumbinary orbits exist for disk material starting at 2-3 
times the binary semi-major axis, depending on eccen- 
tricity (Artymowicz and Lubow 1994). If close binaries 
only impacted a parent disk via inner tidal truncation, we 
would easily detect circumbinary disks extended out to 
~ 100 — 200 au around systems with separations of several 
tens of au, i.e. those near the peak of the binary period distri- 
bution ($2.4). In fact, such disks are uncommon. Czekala 
et al. (2019) compiled a list of known circumbinary disks 
around young stars; among the stars with ages of 1-20 Myr, 
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only 4 of the 30 circumbinary disks are around binaries 
with semi-major axes wider than 10 au. Given the com- 
pleteness and high-sensitivity of ALMA surveys in nearby 
star-forming regions, the lack of circumbinary disks around 
multiples with separations 210-20 au suggests they are 
truly scarce. Anecdotally, at younger ages there a number of 
known circumbinary structures around systems in this sepa- 
ration range ($4.3.2), but it is not clear that the difference in 
disk frequency is statistically significant. If the difference 
is real, it could imply the structures seen in the embedded 
phase are not stable Keplerian disks or that circumbinary 
disks dissipate rapidly during the final phases of binary for- 
mation. As we discuss further in $5.4.2, the reappearance 
of circumbinary debris disks only deepens this puzzle. In 
higher order multiples , it may be that confinement of the 
circumbinary disk due to both the internal binary and an 
external companion slows its viscous evolution and con- 
tributes to a longer lifetime (e.g., HD 98800,0bj]HD 98800 
Ribas et al. 2018; Ronco et al. 2021; TWA 3,0bj]TWA 3 
Czekala et al. 2021). 

Among closer binaries, there may be an anti-correlation 
between the binary mass ratio and the likelihood of retain- 
ing a circumbinary disk. As noted in $2.4.1, in a large sam- 
ple of young stars, Kounkel et al. (2019) find that sources 
with disks are much less likely to be double-lined spectro- 
scopic binaries than those without disks, but that the deficit 
of binaries is recovered when considering single-lined sys- 
tems with RV variability. Thus, perhaps the formation of 
unequal-mass binaries allows or is aided by the presence of 
a circumbinary disk. 


5.3. Orbit-Disk Alignment 


Within a given young binary system, the relative orienta- 
tion of the binary orbital plane and any disks in the system 
provides constraints on binary formation mechanisms (see 
$4.3.4). Even if the birth alignment is substantially altered 
by the time of the Class II phase, disk alignment also pro- 
vides insight into the architecture of future planetary sys- 
tems, since the planets will (at least initially) lie near the 
disk plane. Resolved disk images trace out potential plane- 
tary orbits on the sky and thus provide insights into the ini- 
tial orientation of those orbits, informing our understanding 
of future (mis)alignment. As in the previous sections, we 
separate the discussions of circumbinary and circumstellar 
disks, since their observational constraints differ. 


5.3.1. Alignment of circumbinary disks 


The only systems for which the binary orbital orien- 
tation can be meaningfully constrained on an individual 
(rather than statistical) basis are those that complete a non- 
negligible fraction of their orbits on human timescales, i.e., 
those with periods less than a few hundred years. For low- 
mass stars, this implies semi-major axes less than tens of au 
and thus, as discussed above, any remaining disks (and cer- 
tainly most resolvable disks) tend to be circumbinary rather 
than circumstellar. Despite their rarity, known circumbinary 
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disk systems are valuable since they allow direct compari- 
son of orbital and disk orientations. 

For systems with double-lined spectroscopic orbits and 
circumbinary disks, the binary orbital inclination ?pinary 
can be measured by determining (Mı + M;)sin? Übinary 
from the orbit and the sum of the stellar masses from the ro- 
tation of the disk (Prato et al. 2002; Rosenfeld et al. 2012; 
Czekala et al. 2016, 2017; Prato et al. 2018). For short- 
period binaries, measurements of tpinary generally agree 
with the observed disk inclination taąisk, strongly suggest- 
ing coplanarity of disk and binary orbital planes. Formally, 
tbinary =tdisk does not require coplanarity since OQpinary, 
the position angle of the ascending node of the binary orbit, 
is unobserved and not guaranteed to be the same as the po- 
sition angle of the disk. Nonetheless, Czekala et al. (2019) 
showed observations can constrain the distribution of rela- 
tive orientation angles 0. They found that the data are con- 
sistent with binaries with P <40 days being aligned with 
their disks to within 0 < 5°. Most of these systems also have 
circular orbits. In one case (V4046 Sgrobj]V4046 Sgr), 
coplanarity is confirmed by the detection of eclipse-like 
shadows cast on the disk surface as the stars orbit (D'Orazi 
et al. 2019). 

For longer-period binaries with circumbinary disks, the 
situation is dramatically different; most have eccentric or- 
bits and 0 20? (Czekala et al. 2019). While some re- 
cent discoveries follow these patterns (Czekala et al. 2021; 
Ragusa et al. 2021), notable counterexamples exist, as 
well; WW ChaobjIWW Cha (P=207 days) and V892 
Tauobj]V892 Tau (P — 7.7 yr) are each aligned with their 
circumbinary disks to within ~ 8° (GRAVITY Collaboration 
et al. 2021; Long et al. 2021). Two systems, 99 Herobj]99 
Her (Kennedy et al. 2012) and HD 98800Bobj]HD 98800 
(Kennedy et al. 2019), have disks with a polar orientation to 
the binary orbital plane. 

These two populations agree well with predictions from 
the theory of disk-binary alignment. Short-period systems, 
particularly those that accrete from a circumbinary disk, 
are expected to efficiently align the disk and binary or- 
bital planes on timescales shorter than a typical disk life- 
time (Foucart and Lai 2013). Alignment is faster for short- 
period systems both because it scales with orbital period 
and because these systems circularize their orbits quickly 
(Mathieu 1992); a circular orbit means that the disk in- 
ner edge can be closer to the binary (e.g., Artymowicz and 
Lubow 1994), increasing the alignment torque (Foucart and 
Lai 2014). In contrast, in the presence of sufficient initial 
eccentricity and misalignment, a variant of the Lidov-Kozai 
mechanism (Lidov 1962; Kozai 1962) can cause an outer 
low-mass body (in this case, a disk gas or dust parcel) or- 
biting a massive inner pair to exchange eccentricity for an- 
gular momentum, causing its inclination to librate around 
0 — 90? (Verrier and Evans 2009; Farago and Laskar 2010). 
In a dissipative disk, librations are damped and the disk 
can settle into a polar configuration (Martin and Lubow 
2017). The higher the eccentricity, the lower the initial 
misalignment needed for the disk to reach polar alignment 
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(Martin and Lubow 2018). Some misalignment is required, 
though, providing a constraint on binary formation theo- 
ries; the presence of these polar-aligned disks shows that at 
least some close binaries form with misaligned circumbi- 
nary disks. 

Finally, we note that a circumbinary disk may not be a 
single, planar structure. Theory suggests that the disk can 
be warped or torn into separate rings (Nixon et al. 2013; 
Facchini et al. 2013, 2018), and these effects are observed 
in V892 Tauobj]V892 Tau (Long et al. 2021) and the triple 
system GW Oriobj]GW Ori (Kraus et al. 2020), respec- 
tively. 

In some cases the inner circumstellar disk is mis- 
aligned with the surrounding circumbinary disk obj]HD 
142507(HD 142507; Price et al. 2018; Claudi et al. 2019), 
causing shadowing of the outer disk (Fukagawa et al. 2006; 
Avenhaus et al. 2014). However, it is not clear if all cases 
of inner/outer disk misalignment are caused by stellar com- 
panions since shadows are also present in some systems 
that lack known companions (Benisty et al. 2018; Pinilla 
et al. 2018). 


5.3.2. Alignment of circumstellar disks 


Wider young binary systems, especially those wide 
enough to host individually-resolvable circumstellar disks, 
typically have orbital periods that are too long for their 
orbits to be meaningfully constrained on an individual ba- 
sis. In such systems the orbital inclination is unknown and 
cannot be compared to the disk inclination. However, in 
systems where more than one star has a disk, the orienta- 
tions of the disks can be compared. If the disk inclinations 
are similar, there is a high likelihood that they also lie near 
the binary orbital plane, since the binary orbit has most of 
the system's angular momentum, and thus disk orientations 
tend to be influenced more by the stellar companion than by 
the other star's disk (Bate et al. 2000; Lubow and Ogilvie 
2000; Zanazzi and Lai 2018). 

Itis clear that disks in young binaries are not all aligned, 
as numerous examples exist of individual systems with dif- 
fering disk orientations (Stapelfeldt et al. 1998; Koresko 
1998; Roccatagliata et al. 2011; Jensen and Akeson 2014; 
Williams et al. 2014; Fernández-López et al. 2017). The 
eponymous T Tauri systemobj]T Tauri itself is a clear case 
of misalignment in a triple system (Beck et al. 2020, and ref- 
erences therein). For constraining binary formation models, 
however, systematic studies can inform the extent to which 
disks are aligned, misaligned, or even randomly oriented. 

Some early surveys used polarization of unresolved 
disks to trace disk orientation in binary systems, finding 
that the two stars in most systems have polarization posi- 
tion angles that are more similar than would be expected 
from a random distribution (Monin et al. 1998; Jensen et al. 
2004; Monin et al. 2006). This suggests a tendency toward 
disk alignment, though interpretation of the data is compli- 
cated by the relatively small polarizations in the unresolved 
light and possible contamination by interstellar polariza- 
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tion. Studies with resolved disks are clearly the way for- 
ward, and the advent of ALMA has made it possible to 
resolve disks in many more systems. However, the small 
sizes of disks in binaries (85.1.2) and the need to resolve 
both disks in a given system means that it is still challenging 
to assemble large samples to study this question. 

Jensen et al. (2020) combined new ALMA observations 
with data from the literature to assemble a sample of 8 bi- 
nary systems (projected separations 100—2000 au) in which 
both disks are resolved and disk rotation is detected in CO 
for both components. The latter criterion makes it possi- 
ble to determine whether disks with similar position angles 
and inclinations also rotate in the same direction, allowing 
three-dimensional constraints on the disks’ angular momen- 
tum. Indeed, Jensen et al. (2020) found there is a strong ten- 
dency for both disks in most systems to have similar orien- 
tations, and despite the modest sample size the overall dis- 
tribution is inconsistent with random relative orientations at 
a high level of significance. There is no obvious correlation 
of misalignment angle with projected separation. 


5.3.3. Implications for binary formation models 


One challenge in connecting observed disk alignment to 
binary formation models is the uncertainty in how much 
disk orientations evolve between formation and the 1—2 
Myr ages at which they are observed. Provided there is 
some dissipation in the disk, circumstellar disks should be 
driven toward alignment with the binary orbit (Lubow and 
Ogilvie 2000; Bate 2000; Zanazzi and Lai 2018). This evo- 
lution means that observed Class II samples are likely to 
be more co-aligned on average than the systems were at 
formation. The presence of some systems with misalign- 
ments of many tens of degrees, and few systems that are 
perfectly aligned, thus is consistent with binary formation 
models (e.g. turbulent fragmentation; $3.1) that preferen- 
tially produce systems with disks that are misaligned with 
the binary orbit, at least for systems wider than a few hun- 
dred au. At the same time, theoretical calculations of align- 
ment timescales (e.g., Zanazzi and Lai 2018) predict signif- 
icant tilt evolution in 2 Myr for systems with semi-major 
axes less than ~200 au. Thus, the presence of well-aligned 
systems with separations of order 1000 au and the lack of 
systems with misalignments >90° at any separation sug- 
gest that disk alignment is much faster than currently pre- 
dicted and/or there is some tendency for binary formation 
to favor systems with non-random disk orientations. Rapid 
alignment may be more likely given observations of mis- 
aligned outflows and disks during the protostellar phase 
(84.3.4). Preliminary evidence of a tendency toward align- 
ment between binary companions and the orbits of close- 
in planets around one star (Christian et al. 2022; Dupuy 
et al. 2022) is consistent with the binary-disk alignmentseen 
at young ages, though there may be some difference in 
alignment statistics for close-in terrestrial planets and hot 
Jupiters (Behmard et al. 2022). 

Finally, we note that the tendency toward disk-binary 
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alignment means that most planetary systems in binaries 
will not have a stellar companion with a sufficient relative 
tilt for Kozai-Lidov oscillations to be effective in driving 
high-eccentricity migration, at least not for very long. 


5.4. Constraints From Planets in Multiple Systems 


Theoretical investigations of planet formation and sur- 
vival in binary systems predate their discovery (Harring- 
ton 1977; Szebehely and McKenzie 1981; Innanen et al. 
1997; Holman and Wiegert 1999; Quintana et al. 2002; 
Nelson 2003; Quintana and Lissauer 2006). Early RV sur- 
veys were strongly biased against the detection of planets 
in intermediate and close separation binaries. Nevertheless 
a 2Mjup mass planet was discovered in the y Cephei binary 
in the early 2000s (Hatzes et al. 2003). With a semi-major 
axis of only 2 au in a binary with semi-major axis of 20 au, 
it remains among the tightest binaries with a planet in a cir- 
cumstellar orbit. The first circumbinary planet, Kepler 16b, 
was discovered by the transit method nearly a decade later 
(Doyle et al. 2011). 

Since these first discoveries, there has been a deluge of 
new detections of planets in binaries. Some have derived 
from targeted RV surveys (Eggenberger 2010), coupling 
RV with adaptive optics (Wang et al. 2014, 2015b), Kepler 
and TESS follow-up (Kraus et al. 2016; Wang et al. 2015c; 
Ziegler et al. 2020), and direct imaging (Wang et al. 20152). 
Here, we focus specifically on the statistics and properties 
of planets in multiple systems in so far as they constrain the 
dominant formation mechanism of the host stars. 


5.4.1. 


For circumstellar planets, we have robust data sets 
demonstrating that intermediate-wide separation binaries 
reduce disk masses and truncate disks, thereby substan- 
tially reducing the total mass reservoir available for planet 
formation (see $5.1.1,5.1.2). Planet occurrence rates in 
intermediate separation binaries are correspondingly de- 
pleted compared to single stars. Moe and Kratter (2021) 
estimated that between ^ 10— 200au, planet occurrence 
rates rise from 15% of the single star rate to 100%. It is 
noteworthy that both the disk properties and planet occur- 
rence rates seem to converge to the singleton result around 
200 au. While there are claims in the literature that Hot 
Jupiters reside predominantly in intermediate separation bi- 
naries (Ngo et al. 2016), these claims are not borne out with 
revised statistical analyses (Moe and Kratter 2019). 

The connection between disk properties and planets is 
unsurprising. Binary eccentricity, disk truncation, and over- 
all mass reduction should correlate with a reduced popula- 
tion of observable planets in intermediate separation bina- 
ries. Eccentric binaries especially inhibit planet formation 
by exciting collisional velocities between pebbles and plan- 
etesimals, favoring fragmentation over agglomeration (Sils- 
bee and Rafikov 2015a). Small disks can hinder planet for- 
mation even at disk radii well below the truncation radius; 
outer disks may serve as a vital reservoir for solids that are 
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subject to rapid radial drift (Youdin and Chiang 2004; Na- 
jita and Kenyon 2014). Starving the inner disks of solids 
could also yield planets so small that they are mostly unde- 
tectable (Dupuy et al. 2016). These results imply that plan- 
ets in binaries may have distinct mass and semi-major axis 
distributions. There is only tentative evidence for such a 
shift, data sets only support an overall reduction (Moe and 
Kratter 2021; Hirsch et al. 2021), motivating future dedi- 
cated surveys. 

Extracting constraints on formation channels from inter- 
mediate separation binaries is not straightforward. Disk and 
planet properties both return to their nominal single star val- 
ues at the same separations where we expect disk fragmen- 
tation to be inefficient and turbulent fragmentation to domi- 
nate. It is then tempting to conclude that turbulent fragmen- 
tation does not substantially suppress planet formation. The 
coupled disk and planet trends also suggest that disk mate- 
rial beyond ~ 50 au contributes little to the planetary sys- 
tems that dominate our data sets, i.e., those with orbits < 10 
au. Future Gaia studies of slightly wider planets from e.g., 
astrometric detections, will confirm or refute this claim. 

It is unclear whether a deficit of planets at smaller orbital 
separations indicates a transition between binary formation 
mechanisms. Both turbulent fragmentation with rapid gas- 
driven migration and disk fragmentation could create inhos- 
pitable environments for planets. Thus we believe the cur- 
rent data is insufficient to conclude that either a majority 
of closer separation binaries derive from disk fragmenta- 
tion or that disk fragmentation is inherently less favorable 
to subsequent planet formation. We nevertheless put for- 
ward several hypotheses as to why disk fragmentation may 
hamper planet formation. Binaries formed via disk frag- 
mentation might consume more of the remnant disk mate- 
rial than single stars — both because the instability drives 
rapid accretion prior to fragmentation and because of the 
introduction of a second mass sink. Moreover, as reviewed 
in other PPVII chapters current disk mass measurements in 
single stars point towards the very early growth of solids to 
planetesimal sizes. If gravitational instability prior to frag- 
mentation reduces the survival rate of large grains, these 
systems might never quickly agglomerate dust into planets 
before it is lost to the star. Note that the presence of spi- 
ral arms due to gravitational instability could enhance solid 
growth rates in some cases (Rice et al. 2004; Baehr and Zhu 
2021). Finally, disk substructure around single stars due 
to pressure bumps may be crucial for retaining disk mate- 
rial and growing planets. The maintenance of these delicate 
structures might similarly be inhibited by both an instability 
preceding fragmentation and by the subsequent stirring of 
an embedded, migrating binary. This latter concern would 
be equally valid for binaries that derive from turbulent frag- 
mentation and migrate to intermediate separations. 


5.4.2. Planet occurrence rates in circumbinary disks 


The relative rarity of circumbinary disk detections seems 
at odds with the discoveries from Kepler and TESS (Czekala 
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et al. 2019). The dearth of disks, especially around bina- 
ries with a 10 au is (at least naively) unsurprising, given 
the myriad hurdles to survival posed by close binary for- 
mation models. Dynamical instability, capture, or secular 
oscillations like KL cycles would all be very destructive to 
circumbinary disks (Reipurth and Clarke 2001; Clarke and 
Pringle 1991; Mufioz and Lai 2015; Hamers et al. 2016; 
Martin et al. 2015). Even disk fragmentation and subse- 
quent migration might significantly deplete the circumbi- 
nary mass reservoir. Circumbinary planet statistics are im- 
proving if still limited (Armstrong et al. 2014; Martin 2019; 
Kostov et al. 2020, 2021). The uncertainty is driven not 
only by the ad-hoc nature of most discoveries, e.g., by eye 
rather than via pipeline, though see Martin and Fabrycky 
(2021), but due to unknown inclination distributions. As 
circumbinaries are identified via transit surveys, only those 
well aligned with an eclipsing binary will be easily discov- 
ered. Note that in principal misaligned systems could show 
single transit events. Thus estimates for the occurrence rates 
of planets around close binaries range from 10 — 5096 (Arm- 
strong et al. 2014). 

We cannot easily reconcile these statistics with those 
of circumbinary disks, but ancillary evidence suggests that 
disk detections underestimate the true underlying reservoir 
of mass for forming such planets. Both historical and recent 
surveys of debris disks paint a different picture (Trilling 
et al. 2007; Yelverton et al. 2019). These studies mirror 
the dearth of bright disks in intermediate separation bina- 
ries but find debris disks in close binaries have a frequency 
consistent with that of single stars. One interpretation of the 
dichotomy between protoplanetary disks and debris disks is 
that the lifetime of the bright, gas-rich phase is very short, 
but not so short as to preclude planet formation. 

Whether circumbinary disks provide a net boost or hin- 
drance to planet formation is an area of active inquiry. 
Numerous studies conclude that planetesimal formation is 
more challenging in circumbinary disks due to higher col- 
lisional velocities and nodal precession driven by both the 
disks and binaries (Moriwaki and Nakagawa 2004; Marzari 
et al. 2013; Silsbee and Rafikov 2015b). Bromley and 
Kenyon (2015) have argued that beyond an inner critical 
radius that depends on binary properties, planet formation 
proceeds much like it does around single stars. Formation 
of planetesimals by the streaming instability could amelio- 
rate the bottleneck even at closer separations (Youdin and 
Goodman 2005; Silsbee and Rafikov 2021). If formation 
is inhibited in-situ, inward migration must occur. State-of- 
the-art migration models for more massive planets can often 
reproduce the observed semi-major axis distribution of the 
known planets (Penzlin et al. 2021). This corroborates the 
view that planet scattering and dynamical instability are not 
responsible for the abundance of planets near the stability 
boundary set by the binary orbit (Li et al. 2016; Smullen 
et al. 2016). Though small planets remain undetected as 
of yet, terrestrial planet formation is also expected to pro- 
ceed quickly, at least beyond some inner critical radius 
(Quintana and Lissauer 2006). Recent models by Childs 
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and Martin (2021) found that terrestrial planet formation 
can occur on an accelerated timescale in circumbinary sys- 
tems, although the final masses are typically smaller than 
those in single star systems. This finding is roughly con- 
sistent with the dearth of circumbinary disks. At present all 
known circumbinary planets are above the terrestrial bound- 
ary 6Rearth < R. While smaller planets are harder to find, 
the detection limit is thought to be a factor of two smaller 
(Martin and Fabrycky 2021). Future surveys will clarify 
whether the dearth of small planets is genuine or due to se- 
lection bias. 


6. FUTURE OBSERVATIONS AND OUTLOOK 
6.1. 


While observations and theory have made great progress 
advancing our understanding of the origin of multiplicity 
over the last decade, a variety of open questions remain. 

What portion of multiples form via core versus disk 
fragmentation? Star cluster calculations have demon- 
strated that a single formation mechanism cannot repro- 
duce observed multiplicity statistics. Obtaining statistics 
of forming multiples — catching formation sufficiently early 
that evolution cannot blur the lines between disk and core 
fragmentation — is an observational challenge for the next 
decade. This entails making high-resolution observations 
of Class 0 disks at the earliest possible times and measuring 
their stability, as well as observing cores while collapse is in 
progress. More predictive models of gas-driven migration 
are also necessary to explain the final period distribution, 
especially for the closest binaries. Only a handful of star- 
less cores have confirmed substructure in either dust con- 
tinuum or a dense gas tracer, indicating incipient multiple 
formation (e.g., Kirk et al. 2009; Chen et al. 2010; Naka- 
mura et al. 2012; Lee et al. 2013; Kirk et al. 2017). Future 
surveys of both line and dust emission are needed to iden- 
tify and confirm gas substructures and to ascertain whether 
these peaks represent collapsing regions. 

What are the underlying mechanisms responsible for 
multiple formation in complex numerical simulations? We 
have reviewed the manner in which multi-physics star- 
formation simulations do and do not match observed mul- 
tiplicity statistics. However, the literature often lacks a 
critical analysis of why. "Why did the included physics 
lead to the observed population of multiples? What chan- 
nels of formation dominate on different scales or environ- 
ments? Answering these questions requires tracing the flow 
of mass and angular momentum for each star formed over 
the course of the simulation. While attempts at such analy- 
sis have been made (Smith et al. 2009; Smullen et al. 2020), 
it will require methodologies, like tracer particles, that fol- 
low the time-evolution of gas parcels and capture all impor- 
tant physical processes (e.g., Grudić et al. 2021). Machine 
learning algorithms could be leveraged to identify underly- 
ing evolutionary trends or mechanisms that depart from the 
classical analytic constraints highlighted in $3. 


Summary of Open Questions 
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How do initial conditions, particularly those in extreme 
environments, impact multiplicity? Observations of local 
star-forming regions tentatively find differences in multi- 
plicity properties. Early dynamical evolution, rather than 
innate differences in the gas conditions, appear to drive 
some of these differences (Tobin et al. 2022). Simulations 
suggest that both metallicity and mean magnetic field im- 
pact the primordial multiplicity, however, order of magni- 
tude variations in these conditions may be required to de- 
tect robust statistical variation. Local star-forming regions, 
from Taurus to Orion, exhibit relatively modest physical 
differences. Future studies must push the resolution lim- 
its of statistical surveys or develop new techniques to probe 
multiplicity in more extreme environments. Such observa- 
tions are required to disentangle the influence of nature (ini- 
tial conditions) versus nurture (dynamical evolution). 

Do high-mass stellar multiples form similarly or differ- 
ently from low-mass stars? High-mass protostars (M. = 
8 Mo) are more challenging to observe due to larger 
distances, higher optical depths, and complex chemistry 
(Rosen et al. 2020). Only a few high-mass protobina- 
ries have been observed to date (e.g., Beltrán et al. 2016; 
Beuther et al. 2017b,a; Zhang et al. 2019; Zapata et al. 
2019). Models suggest that high-mass protostellar disks 
are more likely to be gravitationally unstable (Kratter 
and Matzner 2006). Observations, however, are incon- 
clusive (Cesaroni 2005; Chini et al. 2004; Patel et al. 2005; 
Koumpia et al. 2021; Bosco et al. 2019; Motogi et al. 2019). 
Robustly linking O-star multiplicity with high-mass star 
formation requires more high-resolution observations of 
high-mass protobinaries and their disks properties. 

How do the properties and lifetimes of disks in multiple 
systems vary from those in single star systems? — Obser- 
vations demonstrate that disk properties vary with stellar 
mass, age, and environment as well as multiplicity. Larger 
statistical samples are needed to separate these effects and 
determine exactly how disks form and evolve differently in 
multiple systems. To date, star cluster calculations have not 
addressed differences in disk properties in multiple systems. 
Non-ideal magnetic processes and au-scale resolution are 
required for future calculations to adequately follow disk 
formation and evolution. 

What are the separation limits and properties of very 
wide binaries? Previous studies of multiplicity (and 
multiple-system properties) have been limited at wide sep- 
arations by the challenges of confirming whether wide sys- 
tems are bound or only chance projections. With the advent 
of Gaia, it is much more feasible to identify and character- 
ize the widest bound systems (e.g., El-Badry et al. 2021). 
Work to date has only begun to mine the available data. 


6.2. Next Generation Instruments 


A variety of new surveys and observational advances 
over the coming decade hold significant promise for mul- 
tiplicity studies. The next generation Very Large Array 
(ngVLA) is a proposed single-configuration interferometer 
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with ~10 times the resolution and sensitivity of the current 
VLA, which will excel in characterizing close-separation 
protostellar multiplicity in the youngest systems. It is cur- 
rently in the design phase with construction tentatively 
planned to start in 2025. The ngVLA will simultaneously 
offer ~au-scale spatial resolution and high sensitivity, en- 
abling searches for closely separated protostellar multiples 
in nearby star-forming regions at wavelengths of ~7 mm 
and ~3 mm. While ALMA can provide comparable spa- 
tial resolution, the time available for the highest resolution 
studies is quite limited. Meanwhile, the high dust opac- 
ity of protostellar disks at A ~ 1 mm reduces the utility 
of ALMA for characterizing close-separation multiplicity. 
ALMA will soon be able to observe at 8 mm wavelengths, 
though the finest spatial resolution will be less than the 
VLA currently provides. The longer wavelengths of the 
ngVLA will enable it to probe multiplicity within protostel- 
lar disks at wavelengths where the dust is optically thin. 

The higher resolution capabilities of the ngVLA will 
provide insight into multiplicity in more distant massive 
star-forming regions and across more diverse environments. 
The ngVLA sensitivity and resolution in regions a few kpc 
away will be comparable to that of current studies in nearby 
clouds. Since most stars in fact form in high-mass environ- 
ments (e.g., the total number of Class 0 and I protostars 
within 500 pe is ~500 Dunham et al. 2014b), the general 
conclusions that can be drawn from studies of nearby re- 
gions is limited in scope and statistical robustness. 

More distant regions, however, do not currently have 
well-characterized solar-type protostar populations. Accu- 
rate classification is critical to identify YSOs and study 
multiplicity evolution. The James Webb Space Telescope 
(JWST) will provide some of the needed characterization 
for wavelengths less than 24 um. However, photometry 
from the near-infrared to the submillimeter is also required 
and is essential for the youngest sources. Consequently, 
robust measurements of protostellar multiplicity in distant 
star-forming regions will require combining data from next 
generation space-based far-infrared observatories, like the 
proposed Origins Space Telescope, with data from large 
ground-based (sub)millimeter telescopes such as the forth- 
coming Large Millimeter Telescope (LMT) and Fred Young 
Submillimeter Telescope (FYST) or the proposed Atacama 
Large Aperture Submillimeter Telescope (AtLAST). 

Surveys with large ground-based (sub)millimeter tele- 
scopes will enable studies of the environments around 
proto-multiples. The LMT, FYST, and AtLAST telescopes 
will achieve resolutions of <5000 au in dust continuum 
and dust polarization in nearby clouds (< 500 pc). Planned 
surveys with LMT and FYST will resolve large (> 10?) 
core populations, enabling robust comparative studies of 
cores hosting multiple and single systems. Moreover, these 
facilities will better resolve the material surrounding wide 
binaries, and in particular, measure magnetic field proper- 
ties on ~ 10? au scales, constraining the role of magnetic 
fields in fragmentation and dynamical evolution. While an- 
gular resolution will be lower compared to JWST, ALMA, 
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and (ng)VLA making confusion and crowding an issue for 
characterizing multiple systems, they will be able to study 
YSOs in nearby star-forming regions with unprecedented 
angular resolution and observe more distant regions near 
the current spatial resolution. 

Upcoming optical and infrared facilities will probe un- 
charted regions of the parameter space of MS and pre-MS 
binaries (Price- Whelan et al. 2019; Schaefer et al. 2019; Rix 
et al. 2019). JWST, the Roman Space Telescope, GRAV- 
ITY+, and the Giant Magellan Telescope will better char- 
acterize cool, faint BD/late-M binaries and the disk prop- 
erties of T Tauri binaries (Fontanive et al. 2021). Ground- 
based Extremely Large Telescopes will provide the highest 
resolution NIR imaging of star-forming regions and thus 
push multiplicity studies significantly ahead from Class II 
through MS phases. The third full data release of Gaia ex- 
pected in mid-2022 will include multi-epoch spectroscopy, 
photometry, and astrometry of billions of stars, yielding a 
catalog of millions of spectroscopic, eclipsing, astrometric, 
and visual binaries with full orbital solutions (Eyer et al. 
2013, 2015). This catalog will likely reveal hidden trends 
in multiplicity statistics, such as the relationship between 
stellar mass and the distributions of triple star periods, mass 
ratios, and orientations. Finally, the Large Synoptic Survey 
Telescope at the Vera Rubin Observatory will discover tens 
of millions of eclipsing binaries (Prsa et al. 2011; Geller 
et al. 2021a), providing insight into how close binary prop- 
erties vary with stellar properties and environment. 


6.3. Challenges for Theoretical Models 


While current state-of-the-art simulations reproduce a 
number of key observables, such as the strong dependence 
of multiplicity on stellar mass ($3.4), most do not repro- 
duce the observed separation distribution or its evolution 
over time. In addition, models must also reproduce the ob- 
served triple fraction and companion frequency as shown in 
Fig. 1. Published comparisons often benchmark numerical 
results using MS field star population statistics, which are 
inappropriate for simulated sources < 0.5 Myr years old. 
It is critical to make equivalent comparisons by account- 
ing for observational constraints, including sensitivity lim- 
its and projected separation (rather than 3D separation or 
binding energy). Finally, public release of numerical meth- 
ods and simulation data will increase transparency and en- 
courage broader community participation. 

To obtain robust stellar statistics, future numerical stud- 
ies must capture a dynamic range approaching seven orders 
of magnitude in spatial scale, while also including all im- 
portant physical processes. This requires exploiting com- 
putational advances in clock speeds and system architec- 
tures and developing novel, efficient methods to model star 
cluster formation. Future star-cluster calculations must also 
incorporate non-ideal magnetic effects, which are important 
for modelling disk evolution. Including these processes, to- 
gether with sub-au resolution, may mitigate current discrep- 
ancies with observation and point the way to a fully realized 
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theory of multiple formation. 


7. CONCLUSIONS 


We have reviewed the formation and evolution of stel- 
lar multiplicity from protostellar birth through the main se- 
quence lives of stars. Observations and theory agree that 
most stars are born in small-order multiple systems. Multi- 
plicity arises primarily from the fragmentation of filaments, 
dense cores and massive accretion disks. Young stellar sys- 
tems evolve through dynamical interactions such that mul- 
tiplicity declines with time. System statistics, such as the 
fraction of multiple systems, frequency of higher order mul- 
tiples and median separation, are strong functions of pri- 
mary stellar mass. These also vary with environment den- 
sity and metallicity. 

While angular momentum misalignment is common dur- 
ing the embedded phase, as indicated by protostellar out- 
flow and disk orientations, older systems more often exhibit 
orbit-spin alignment and have aligned disks. Protoplane- 
tary disks in multiple systems are on average less massive 
and more compact than those around single stars, while cir- 
cumbinary protoplanetary disks are rare. These trends sug- 
gest planet formation and planetary architectures differ in 
multiple systems. Additional high-resolution surveys with 
existing facilities, such as ALMA and Gaia, as well as data 
from future facilities such as the JVLA, extended ALMA 
and JWST, promise significant strides towards addressing 
the remaining mysteries of multiple star systems. 
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